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INTERPLANETARY  MAGNETIC  FIELD  AND  TROPOSPHERIC  CIRCULATION 


John  M.  Wilcox ,  Philip  H.  Scherrer  and  J.  Todd  Hoeksema 
Institute  £or  Plasma  Research,  Stanford  University 
Stanford,  California 


ABSTRACT 

The  relation  between  interplanetary  magnetic  sector  boundary 
crossings  and  areas  of  high  vorticity  in  the  troposphere  that  was 
reported  during  1963-1973  cannot  be  investigated  in  the  years 
after  1973  because  of  changes  in  the  processing  of  the  500  mb 
height  grids  prepared  by  the  National  Meteorological  Center.  In 
particular,  we  cannot  say  that  the  effect  disappeared. 

The  same  applies  to  vorticity  computed  from  NMC  winds  grids. 

The  Limited  Area  Fine  Mesh  grid  has  a  large  noise  in  com¬ 
puted  vorticity  after  December  3,  1974.  Therefore  the  interest¬ 
ing  analysis  of  Larsen  and  Kelley  cannot  be  extended.  They  had 
found  that  forecasts  of  Vorticity  Area  Index  were  significantly 
poorer  after  a  sector  boundary. 


1.  INTRODUCTION 


We  have  reported  (Wilcox,  et  _al.  1973)  that  an  index  that 
measures  the  area  of  intense  circulation  in  the  wintertime  north¬ 
ern  hemisphere  has  a  minimum  value  about  one  day  after  the  solar 
wind  has  carried  an  interplanetary  sector  boundary  past  the 
earth.  This  Vorticity  Area  Index  (VAI)  was  defined  by  Roberts 
and  Olson  (1973)  as  the  area  within  the  vorticity  contours  of  20 
x  10™5  s-1  plus  the  area  within  contours  of  24  of  the  same  units. 
The  second  contour  level  was  added  to  increase  the  influence  of 
the  regions  of  most  intense  circulation.  The  absolute  vorticity 
was  computed  from  the  National  Meteorological  Center  500  mb 
height  grids  for  the  northern  hemisphere  north  of  20  degrees  N, 
using  a  simple  five-point  height  difference  computation.  These 
grids  were  archived  by  the  National  Center  for  Atmospheric 
Research  in  Boulder. 

This  minimum  in  the  VAI  was  discovered  following  54  sector 
boundaries  during  the  winter  months  November-March  of  1964-1970 
(Wilcox  et  _al.  1974).  A  later  analysis  (Wilcox  et  al.  1976) 
using  131  boundaries  in  the  winters  1963-1973  found  a  similar 
minimum,  as  shown  in  Figure  1.  Hines  and  Halevy  (1975,  1977) 
extensively  analyzed  our  reports  and  made  the  following  state¬ 
ment:  "In  reaching  our  initial  decision  to  accept  the  alleged 
solar  signal  as  a  physically  meaningful  signal  and  in  holding  to 
that  decision  we  are  not  asserting  that  physical  reality  has  been 
established  beyond  all  question.  We  assert  only  that  the  demands 


of  the  scientific  method  have  been  met  within  what  we  consider  to 
be  reasonable  limits  and  that  the  time  has  come,  for  us  at  least, 
to  accept  the  physical  reality  of  the  solar  signal  as  what  might 
be  called  the  most  appropriate,  if  not  the  only  appropriate, 

working  hypothesis  of  the  moment.  As  before,  we  invite  our 

readers  to  share  this  view."  (Hines  and  Halevy,  1977). 

An  example  of  the  amount  of  organization  that  existed  during 
the  winters  1963-1973  is  shown  in  Figure  2  (Wilcox,  1979).  The 
northern  hemisphere  was  divided  into  seven  zones  with  approxi¬ 
mately  the  same  area  in  each  zone.  Within  each  zone  an  analysis 
similar  to  that  shown  in  Figure  1  was  computed  except  that  the 
VAI  was  computed  as  the  area  within  vorticity  contours  of  a  sin¬ 
gle  value,  which  was  in  successive  analyses  20,  22,  24... 30  x 
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10  s  .  The  fractional  effect  shown  as  the  ordinate  in  Figure 
2  is  simply  the  depth  of  the  minimum  in  Figure  1  divided  by  the 

average  value  during  adjacent  days.  The  zones  are  largely 

independent  since  an  individual  region  of  high  vorticity  would 
usually  be  in  one  or  at  most  two  zones.  The  minimum  was  found  in 
each  of  these  zones,  and  within  each  zone  the  areas  with  more 
intense  circulation  (represented  by  larger  values  of  vorticity) 
had  a  deeper  minimum  after  the  sector  boundary  transit.  This 
amount  of  organization  appears  unlikely  to  be  associated  with  a 
statistical  excursion. 


Temporal  Changes  in  VAI  Computed  from  NMC  Data 


Williams  and  Gerety  (1978)  claimed  that  the  minimum  dis¬ 
cussed  above  diminished  or  disappeared  during  the  winters  1974- 
1976  that  followed  the  original  analysis.  In  response  we  cited  a 
discovery  by  Hines  and  Halevy  (1977).  They  defined  for  each 
boundary  an  "Excursion"  as  being  the  range  in  VM  (maximum  value 
minus  minimum  value)  during  a  12  day  window  centered  on  the  time 
of  the  boundary.  As  shown  in  Figure  3  (Wilcox  and  Scherrer, 
1981),  boundary  crossings  that  occurred  during  times  of  larger 
Excursions  were  followed  by  a  minimum  whose  depth  did  not  change 
with  time  from  1965  through  1977.  Boundary  crossings  that 
occurred  during  times  of  smaller  Excursions  showed  a  null  effect 
during  most  of  the  years  shown  in  Figure  3. 

Why  then  did  Williams  and  Gerety  (1978)  report  that  the 
effect  had  diminished?  The  explanation  is  found  in  Figure  4 
(Wilcox  and  Scherrer,  1981),  which  shows  for  each  winter  between 
1963  and  1978  the  number  of  boundary  crossings  that  occurred  dur¬ 
ing  times  of  larger  Excursions  and  during  of  times  of  smaller 
Excursions.  During  the  interval  1963-1973  that  was  discussed  by 
Wilcox  £t  aJL.  (1976)  the  number  of  boundaries  with  larger  Excur¬ 
sions  was  approximately  equal  to  the  number  of  boundaries  with 
smaller  Excursions.  After  1974  this  situation  changed  dramati¬ 
cally,  and  in  1977  there  were  no  boundaries  with  larger  Excursion 
and  in  1978  only  one  such  boundary. 


Now  that  a  few  more  years  have  gone  by  it  is  possible  to  get 
a  broader  perspective  on  the  change  with  time  of  the  VAI  as  com¬ 
puted  from  the  National  Meteorological  Center  500  mb  height  grid, 
figure  5  shows  a  plot  of  VAI  from  1964  through  1981.  In  the 
first  years  the  average  value  of  the  VAI  is  larger,  and  a  very 
clear  seasonal  variation  is  evident,  both  in  the  maximum  values 
of  the  index  and  in  the  minimum  values  of  the  index.  In  the 
latter  half  of  Figure  5  a  steady  decrease  in  average  value  of  the 
VAI  is  evident,  and  the  seasonal  variation  that  was  so  clear  in 
the  early  years  is  much  less  evident.  In  fact  during  the  last 
years  shown  in  Figure  5  the  seasonal  variation  in  the  minimum  of 
the  VAI  is  almost  absent.  Our  discussions  with  the  National 
Meteorological  Center  have  revealed  that  in  recent  years  the  500 
mb  height  grid  has  been  subjected  to  increasing  amounts  of  spa¬ 
tial  smoothing  in  order  to  remove  noise.  "There  is  ample  reason 
to  believe  that,  the  data  sat  (500  mb  height  and  winds)  may  well 
be  inhomogeneous  for  the  VAI  computation"  (Cerrity,  1982). 

We  believe  that  this  increased  smoothing  during  the  data 
processing  has  removed  completely  the  intervals  of  large  Excur¬ 
sion  that  were  shown  in  Figure  3  to  be  associated  with  the 
effect.  I t  is  therefore  not  possible  to  make  any  statement  about 
whether  or  not  the  effect  reported  by  Wilcox  et  al .  (1973)  con¬ 
tinued  in  the  years  after  the  original  analysis. 


Figure  6a  shows  an  analysis  similar  to  Figure  1  for  the 
winters  1963-1973,  and  Figure  6b  shows  the  same  analysis  for  the 


winters  1979-1982.  Figure  6b  shows  a  null  effect  in  these  years, 
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but  the  average  value  of  the  VAI  is  only  30  x  10  km  as  compared 

5  2 

with  an  average  value  of  50  x  10  km  in  Figure  6a.  It  is 
exceedingly  unlikely  that  the  average  area  of  high-vorticity 
regions  is  now  only  3/5  as  large  as  it  was  ten  years  ago. 

Since  it  is  also  possible  to  compute  the  VAI  directly  from 
winds,  we  obtained  the  National  Meteorological  Center  500  mb 
winds  grids  in  digital  form  as  archived  by  the  National  Center 
for  Atmospheric  Research  in  the  hope  that  perhaps  less  smoothing 
had  been  applied  to  this  data  set.  Figure  7  shows  the  VAI  com¬ 
puted  from  winds  as  a  function  of  time  from  1965  through  1981. 
We  see  that  the  average  VAI  in  1972  was  about  half  the  average  in 
1971  and  before,  and  that  the  years  after  1972  average  less  than 
half  as  much  as  the  years  before  1972.  Again  the  seasonal  varia¬ 
tion  in  both  maxima  and  minima  is  very  apparent  before  1972  and 
much  less  prominent  after.  This  winds  grid  data  set  is  thus  also 
inhomogeneous  for  VAI  computations. 


3.  Noise  in  the  VAI  Computed  from  LFM  Data 


Larsen  and  Kelley  (1977)  supported  the  result  shown  in  Fig¬ 
ure  1  in  an  analysis  using  the  Limited  Area  Fine  Mesh  (LFM)  Grid 
that  includes  most  of  the  North  American  continent.  They  further 
reported  that  on  the  two  days  following  a  sector  boundary  cross¬ 
ing  the  accuracy  of  forecasts  of  the  VAI  was  considerably  less 
than  average.  Larsen  and  Kelley  analyzed  the  winters  1972  to 


In  the  hope  of  extending  this  interesting  investigation  into 
later  years  we  made  use  of  the  LFM  grid  as  archived  in  digital 
form  by  the  National  Canter  for  Atmospheric  Research.  We 
discovered  that  a  change  in  the  method  of  processing  the  LFM  data 
on  December  2,  1974  introduced  noise  into  all  subsequent  computa¬ 
tions  of  the  VAI  to  such  an  extent  that  it  is  impossible  to 
extend  the  investigations  of  Larsen  and  Kelly.  The  nature  of 
this  noise  is  illustrated  in  the  following  figures:  Figure  8a 
shows  500  mb  height  contours  for  December  1,  1974,  the  last  day 
that  was  free  of  the  noise.  Figure  8b  shows  the  corresponding 
vorticity  contours.  Figure  9a  shows  500  mb  height  contours  for 
December  3,  1974.  Note  the  small  zig-zags  in  the  encircled 
height  contour.  These  lead  to  the  diamond-like  noise  contours 
when  vorticity  is  computed,  as  shown  in  Figure  9b.  Each  LFM  grid 
from  December  3  to  the  present  time  shows  similar  noise  charac¬ 
teristics.  Figure  10  is  the  average  vorticity  computed  for  the 
month  of  November  1975.  The  largest  amplitude  noise  is  in  the 
approximate  geographical  locations  of  the  Rocky  Mountains  and  the 
Sierra  Nevada  Mountains. 

We  have  been  informed  by  the  National  Meteorological  Center 
that  the  onset  of  this  noise  on  and  after  December  3,  1974  can  be 
identified  with  a  change  in  their  processing  procedure  on  that 
date . 
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FIGURE  CAPTIONS 


Figure  1. 


Figure  2. 


Figure  3. 


Figure  4. 


Superposed  epoch  analysis  of  the  500  mb  vorticity 
area  index  (the  area  of  low  pressure  troughs  in  the 
Northern  Hemisphere)  about  times  when  solar  magnetic 
sector  boundaries  were  carried  past  the  earth  by  the 
solar  wind.  The  VAI  is  influenced  by  the  large-scale 
solar  sector  structure,  not  by  the  boundary  passages 
as  such.  Figure  la  uses  50  of  the  54  boundaries  used 
in  the  original  work.  Figure  lb  uses  31  new  boundary 
passages  not  included  in  the  original  analysis,  and 
Figure  lc  is  a  subset  of  Figure  lb  in  which  the 
time  of  46  boundary  passages  were  determined  from 
spacecraft  observations  (Wilcox  et  al,  1976). 

The  fractional  depth  of  the  sector- related  minimum  in 
vorticity  area  index  as  a  function  of  the  vorticity 
value  used  in  computing  the  index  (Wilcox,  1979). 

The  size  of  the  Sun-Weather  effect  for  the  groups  of 
boundary  transits  having  larger  Excursions  (open 
circles)  and  for  the  groups  of  boundary  transits  having 
smaller  Excursions  (filled  circles).  The  total  length 
of  the  error  bar  is  twice  the  standard  error  of  the 
mean  (Wilcox  and  Scherrer,  1981). 

The  number  of  boundary  transits  in  each  year  for  which 
the  Excursions  are  in  the  larger  group  (open  circles) 


0- 


and  in  the  smaller  group  (filled  circles).  Note  that 
the  intensity  of  tropospheric  circulation  as  measured 
by  the  VAI  declined  after  1974  and  decreased  very  much 
after  1976  (Wilcox  and  Scherrer,  1981). 

Figure  5.  The  VAI  computed  from  the  National  Meteorological 
Center  height  grid  at  500  mb  as  a  function  of  time 
from  1964  through  1981. 

Figure  6a.  Similar  to  Figure  1,  showing  a  superposed  epoch 

analysis  of  the  500  mb  VAI  about  sector  boundary 

transits  in  the  winters  1964-1973.  Note  that  the 
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average  value  of  the  VAI  is  about  50  x  10  km  . 

Figure  6b.  The  same  as  Figure  6a,  except  for  the  winters 

1979-1982.  The  minimum  in  VAI  after  sector  boundary 

transit  that  was  so  persistent  in  the  original 

analysis  is  not  present,  but  note  that  the  average 
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value  of  the  VAI  is  now  only  about  30  x  10  km  . 

Figure  7.  The  VAI  computed  from  the  National  Meteorological 

Center  winds  grid  at  500  mb  from  1965  through  1981. 

Figure  8a.  Height  contours  at  500  mb  in  the  Limited  Area  Fine 
Mesh  Grid  of  the  National  Meteorological  Center  for 
December  1,  1974  at  0  Z. 

Figure  8b.  Contours  of  vorticity  calculated  from  Figure  8a. 


Figure  9a 


The  same  as  Figure  3a,  but  for  December  3,  1974  0  Z 


Note  the  zig-zags  in  the  encircled  contour. 

Contours  of  vorticity  calculated  from  Figure  9a. 

Note  the  diamond- shaped  noise  contours  that  are  most 
prominent  on  the  left  hand  half  of  the  figure. 

The  same  as  Figure  9b,  but  averaged  over  the  month  of 
November  1975  to  bring  out  the  noise  contours  more 
clearly. 
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Fig.  10 
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ABSTRACT 

A  numerical  algorithm  using  implicit  time-differencing  is  applied  to  the  solar  wind  equations 
allowing,  for  the  first  time,  solutions  including  thermal  conduction  to  be  found  for  time-dependent 
flow  traversing  the  subsonic  to  supersonic  velocity  transition  region.  Sample  solutions  are  shown  that 
demonstrate  the  distinctive  differences  introduced  by  including  thermal  conduction,  in  comparison  to 
the  commonly  available  solutions  assuming  polytropic  flow.  Also,  it  is  found  that  steady  solutions  are 
produced  at  least  as  quickly  using  a  time-dependent  relaxation  to  the  steady  state  as  when  solving  the 
steady-state  equations. 

Subject  headings:  hydromagnetics  —  Sun:  corona  —  Sun:  solar  wind 

I.  INTRODUCTION 

Numerical  solutions  of  the  time-dependent  solar  wind  equations  have  been  appearing  for  several  years,  usually  in  the 
context  of  modeling  a  transient  effect  in  the  solar  atmosphere  or  interplanetary  medium.  In  order  to  avoid  problems 
associated  with  treatment  of  boundary  conditions  near  the  Sun,  in  a  subsonic  flow  regime,  most  of  these  solutions  have 
been  confined  to  regions  where  the  flow  is  strictly  supersonic  (see,  e.g.,  Hundhausen  and  Gentry  1969;  Steinolfson, 
Dryer,  and  Nakagawa  1975;  Metzler  etal.  1979;  Pizzo  1980).  These  solutions  have  nevertheless  dealt  with  several 
important  physical  processes  in  the  interplanetary  medium,  such  as  magnetic  field  effects  (Steinolfson,  Dryer,  and 
Nakagawa  1975)  and  multiple  species  interactions  with  thermal  conduction  (Metzler  etal.  1979).  In  contrast,  other 
solutions  have  specifically  addressed  transient  processes  in  the  corona,  dealing  with  the  boundary  condition  problem  in 
two  different  ways.  One  approach  has  been  to  take  an  initially  hydrostatic  atmosphere  with  nonreflective  boundary 
conditions  at  the  outer  boundary  (see,  e.g.,  Han,  Wu,  and  Nakagawa  1981  for  a  very  comprehensive  model  using  this 
approach).  In  a  second  approach,  the  boundary  conditions  at  the  base  in  the  presence  of  an  expanding  corona  have 
been  treated  in  a  self-consistent  manner  (Nakagawa  and  Steinolfson  1976;  Steinolfson  and  Nakagawa  1976).  This 
normally  leads  to  more  mathematical  complexity,  but  has  been  found  to  be  necessary  in  view  of  the  fact  that  the 
corona  cannot  be  approximated  as  hydrostatic  above  a  fraction  of  a  solar  radius  beyond  the  photosphere. 

An  important  limitation  in  the  analyses  of  Nakagawa  and  Steinolfson  (1976)  and  Steinolfson  and  Nakagawa  (1976) 
is  their  assumption  of  adiabatic  flow.  Here,  an  analysis  is  described  which  is  specifically  designed  to  overcome  this 
limitation  and  to  improve  the  efficiency  with  which  both  steady  and  transient  solutions  to  the  solar  wind  equations  can 
be  found.  Specifically,  an  implicit  (see,  e.g.,  Richtmyer  and  Morton  1967  for  a  definition  of  the  terms  explicit  and 
implicit )  finite  difference  scheme  is  developed  for  the  numerical  solution  of  the  solar  wind  equations  including  thermal 
conduction,  for  cases  where  the  flow  is  subsonic  near  the  Sun  and  supersonic  far  from  the  Sun.  It  is  the  intent  in  this 
paper  to  describe  the  analysis  and  numerical  algorithm  in  some  detail  and  to  present  results  comparing  transient  effects 
in  the  corona  for  polytropic  flow  with  those  for  thermally  conductive  flow. 

To  the  best  of  this  author’s  knowledge,  implicit  differencing  schemes  have  been  applied  to  the  solar  wind  equations 
in  only  two  published  instances.  In  the  first  case  Metzler  etal.  (1979)  used  this  approach  to  solve  the  time-dependent 
two-fluid  equations  in  the  supersonic  flow  regime.  In  the  second,  Han,  Wu,  and  Nakagawa  (1981)  analyzed  coronal 
transients  on  a  hydrostatic  initial  state,  in  two  dimensions,  with  a  magnetic  field.  Metzler  etal.  used  implicit 
differencing  in  order  to  be  able  to  include  thermal  diffusion,  whereas  Han  etal.  although  initially  including  thermal 
diffusion,  only  showed  results  for  a  transient  in  a  polytropic  gas.  The  reason  Han  etal.  turned  to  implicit  time 
differencing  was  that  the  limitation  on  step  size  due  to  the  Courant  condition  (Richtmyer  and  Morton  1967)  that 
applies  to  all  explicit  schemes  was  causing  an  impossibly  small  time  step  to  be  used  for  only  moderately  large  magnetic 
field  strengths. 

In  general,  the  solar  wind  equations  are  prime  candidates  for  solution  with  implicit  timc-diffeiencing  because  they 
are  extremely  “stiff’— that  is,  the  propagation  speed  of  a  thermal  pulse  is  much  larger  than  the  sound  or  Alfven 
speeds.  This  can  be  demonstrated  by  comparing  the  characteristic  time  scales  for  propagation  of  a  sound  wave  ( r, )  and 
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a  thermal  wave  (r,.).  These  are  (Craig  and  McClymont  1976): 

r,*  8.3X10 -SL/Ti/2,  rt.«:  3  X  10“  wnL2/Ts/1 , 

where  L,n,  and  T  are  the  length  scale  (cm),  electron  number  density  (cm-3),  and  temperature.  In  the  corona,  L  ~  10" 
cm.  T  ~ 106  K,  and  n  ~  10s  cm-3.  For  these  numbers,  the  ratio  (t,/tc.)  is  about  25,  a  large  number  for  these  purposes. 
This  is  the  reason  for  the  relative  importance  of  thermal  conduction  in  the  present  analysis  in  comparison  to  analyses 
describing  transient  motion  near  the  transition  region  (Craig  and  McClymont  1976;  Wu  et  al.  1981)  where  the  density 
is  much  larger.  Implicit  schemes  are  not  subject  to  the  Courant  numerical  stability  condition  set  by  the  largest 
characteristic  speed,  and  so  are  stable  for  arbitrarily  large  time  steps.  The  penalty  for  this  stability  is  that  implicit 
schemes  are  algebraically  complex  and  hence  have  a  large  arithmetic  operation  count. 

The  above-mentioned  applications  of  implicit  time-differencing  to  the  solar  wind  equations  used  classical  algorithms 
developed  early  in  the  history  of  implicit  schemes,  and  suffer  as  a  consequence.  Han,  Wu,  and  Nakagawa  (1981)  utilize 
the  implicit  continuous-fluid  Eulerian  (ICE)  method  (see,  e.g.,  Harlow  1973).  The  ICE  method,  although  successfully 
applied  to  a  large  class  of  problems,  is  relatively  cumbersome  due  to  its  use  of  a  nonsimultaneous  solution  of  the 
equations  on  a  split  grid.  Metzler  et  al.  (1980)  solve  the  equations  in  a  Lagrangian  frame  of  reference,  treating  the 
energy  equation  implicitly  but  the  momentum  equations  explicitly.  Additional  complexity  results  from  an  iterative 
technique  utilized  to  account  for  the  nonlinearity  introduced  by  the  transport  coefficients.  Recently,  implicit  methods 
have  gained  new  stature  from  the  development  of  concise,  noniterative,  one-dimensional  and  alternating  direction 
implicit  (ADI)  multidimensional  schemes  by  Lindemuth  and  Kileen  (1973),  Briley  and  McDonald  (1975),  and  Beam 
and  Warming  (1978),  among  others.  The  paper  by  Beam  and  Warming  (1978)  will  be  used  extensively  in  the  discussion 
that  follows,  and  so  will  simply  be  referred  to  as  BW. 

The  BW  algorithm,  as  applied  here,  is  noniterative  and  retains  as  nearly  a  conservation  law  form  as  is  possible  for 
the  solar  wind  equations.  This  is  essential,  together  with  the  inclusion  of  dissipation,  for  the  proper  treatment  of 
embedded  shock  waves  (“shock  capturing”).  The  temporal  difference  approximation  is  generalized  to  retain  a  variety 
of  schemes  including  a  three-level  scheme  requiring  only  two  levels  of  data  storage.  The  development  makes  extensive 
use  of  the  BW  “delta”  form  (increments  of  the  conserved  variable  and  flux  vectors)  to  achieve  numerical  efficiency  and 
retain  the  property  of  a  steady  state  (if  one  exists)  independent  of  the  time  step. 

In  the  following  section,  the  physics  of  the  model  and  the  equations  of  motion  are  stated  and  discussed.  The  implicit 
algorithm  is  developed  in  §  III,  and  a  description  is  given  of  how  artificial  dissipation  should  be  added  to  the  algorithm 
in  those  cases  for  which  it  is  found  necessary  (e.g.,  shock  propagation)  in  §  IV.  Section  V  contains  a  discussion  of  the 
boundary  condition  treatment  for  steady  and  time-dependent  flows,  and  examples  are  described  in  §  VI.  Additions  to 
the  equations  and  modifications  of  the  algorithm  to  address  important  questions  in  solar  and  stellar  wind  studies  are 
discussed  in  §  VII,  and  the  results  are  summarized  in  §  VIII. 


II.  MODEL  DESCRIPTION 

The  model  equations  are  for  a  radial  (spherically  symmetric)  flow,  one-fluid  description  of  the  proton-electron  solar 
wind  with  no  momentum  or  energy  addition,  and  no  “superradial  divergence.”  This  set  of  equations  is  sufficiently 
comprehensive  to  produce  interesting  new  results  in  modeling  transient  processes.  The  mass,  momentum,  and  energy 
conservation  equations  for  this  ensemble  of  assumptions  and  model  descriptions  are  then  usually  written: 

|+^(rW,=0.  a) 

(  9ur  dvr  \  _  9 p  CM  ,,, 

p(-97  +  u'07)-_97'p^~’  (2) 

Ji  +  72Yrirvr(e  +  p)]+-rr(rg)^Pvr—=0.  (3) 

where  p  is  the  density,  vr  is  the  radial  velocity,/?  is  the  pressure,  T is  the  temperature,  q  is  the  heat  flux  density,  e  is  the 

total  (kinetic  plus  internal)  energy  per  unit  volume,  and  M  is  the  mass  of  the  Sun.  The  density  is  related  to  the  number 
density  of  protons  and  electrons  by 

p-n(m„  +  mr). 


(4) 
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and  the  perfect  gas  law  is  used: 


p  -  2nkT=2p  j ) T  =  2pRT , 


R  *  8.317X  107  (cgs).  (5) 

As  the  intent  here  is  to  compare  polytropic  flow  with  the  consequences  of  the  presence  of  strictly  classical  conduction, 
the  form  chosen  for  q ,  the  heat  flux  density,  is  the  collision  dominated  model  described  by  Spitzer  ( 1 962): 

q=-KT^fr,  (6) 

1.99X10-5  ,  , 


A  *=  1.246X  104(7’3/«),/2.  (8) 

Equations  (6)— (8)  differ  from  analogous  choices  in  some  other  models  (e.g.,  Holzer  and  Leer  1981;  Whang  and  Chang 
1965)  in  including  the  spatial  variation  of  the  logarithmic  term  in  (8).  This  is  usually  a  small  effect,  but  it  is  easily 
included  in  the  analysis  and  so  shall  be  retained.  Because  this  collisional  model  for  q  is  valid  only  near  the  base  of  the 
corona,  the  computational  examples  in  section  6  will  not  be  carried  beyond  about  20  solar  radii. 

Putting  equations  ( 1)— (3)  into  near  “conservation  law”  form  gives: 

^  +  ^=0  (9) 

3/  3 r  '  { ’ 


£+£(£+,v)-2„-^=0. 


b  =  pr2 


mr  =  pvrr-. 


are  the  conserved  variables  related  to  the  mass,  momentum,  and  energy  densities,  respectively.  The  equation  of  state  is 


.  _  Y  —  1 


e“2T 


where  y  is  the  ratio  of  specific  heats  (or  the  polytropic  index  in  case  thermal  conduction  is  not  included).  The 
temperature  is  given  by 


7'==2S^(T‘',('i’“5m') 


Equations  (9)— ( 11)  can  also  be  written  in  the  vector  form: 


+wi+ll=0. 
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where  the  vectors  U,  F,  V ,  and  H  can  be  deduced  from  (9)-(ll)  by  inspection  (see  Appendix,  [A1]-[A4])  and 
Ur  =  dU/dr.  This  vector  equation,  together  with  boundary  conditions  to  be  discussed  later,  constitutes  the  complete 
description  of  the  problem. 

HI.  algorithm  for  solution  by  implicit  time-differencing 

In  (15),  U  is  the  vector  of  conserved  variables,  F  and  V  are  flux  vectors,  and  H  is  the  inhomogeneous  part  of  the 
equation.  The  treatment  of  this  equation  in  this  section,  along  with  the  development  of  the  implicit  algorithm,  will 
closely  follow  the  discussion  in  BW,  to  which  the  reader  is  referred  for  more  details. 

For  advancing  the  solution  in  time,  BW  suggest  a  general  single-step  temporal  scheme  that  can  be  written  in  the 
form: 


v  =  v[  £  4l,) 


and  AU"  —  U"+'  -U". 


Here,  except  for  the  first  time  step,  X  =  1  and  £  =  { ,  which  produces  a  three-level  second  order  accurate  scheme.  At  the 
first  time  step,  because  only  one  level  of  information  will  be  available,  X  =  ^  and  £  =  0. 

If  (15)  is  solved  for  9l//3i  and  the  resulting  expression  for  the  temporal  derivative  is  inserted  in  (16),  then  it  can  be 
shown  that 


At/"  =  - 


3A  F-  dAV"  1  At 
~Tr - J  +  T+i 


3 Fn  dP" 


+  j^Al/"-'  +  o[(\-}-i)At2  +  Ar3], 


F"  —  FI  2  Ar,  ,  A Fn  =  Fn+'-F\ 


with  similar  definitions  for  AF"  and  AH".  Henceforth,  having  served  its  purpose,  the  order  symbol  0[(X  —  j  ~  £)A t2 
+  At3]  will  be  dropped.  A  problem  arises  because  the  vector  increments  A F",  AV"  and  AH"  are  nonlinear  functions  of 
the  conserved  variables,  V.  A  linear  equation  with  the  same  temporal  accuracy  as  (17)  can  be  obtained  if  a  Taylor 
series  expansion  is  used  (Briley  and  McDonald  1973;  Beam  and  Warming  1976),  such  that 


/•»+i  =  *•<•  + j"(t/"+' -£/»)  + o[A/2] 


AF"  =A"AU"  +  0[At2],  (19a) 

where  A  is,  in  tensor  notation,  the  Jacobian  matrix  3F/31/  (Appendix,  eq.  (A5]).  Significantly,  this  technique  is 
described  as  the  “most  robust  strategy  for  linearization”  by  Dwyer  (1981).  The  terms  in  V  and  H  can  be  treated 
similarly  to  give 

AV"  =  {C-Pr)AU"  +  ^-{PAU)"  +  0[At2],  (19b) 

dr 

AH"  =  EAU"  +  0[At2],  (19c) 

where  Cis  (in  tensor  notation)  the  Jacobian  matrix  dV/dU,  E  is  the  Jacobian  matrix  3/// 31/,  P  is  the  Jacobian  matrix 
W/W„  Pr  =  dP/dr,  and  Ur  =  dU/dr.  These  matrices  are  shown  in  detail  in  the  Appendix  (eqs.  [A6]-[A9]). 
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Introducing  the  approximations  (19)  into  (17)  gives 


where  here,  and  throughout  the  remainder  of  the  paper,  notation  of  the  form 


denotes 


j-r(A  +  C-P,)\w' 


\(A+C-Pr)  At/"],  etc. 


In  (20),  the  right  hand  side  of  the  equation  can  be  recognized  as  the  steady-state  portion  of  equations  (l)-(3)  plus  a 
term  depending  on  A Un~ '.  These  parts  of  (20)  only  depend  on  information  from  the  present  and  previous  time-steps, 
and  so  are  called  the  “explicit”  portion  of  the  equation.  Conversely,  the  left  hand  side  of  (20)  is  all  multiplied  by  A Ua, 
thereby  depending  on  information  at  the  advance  time-step,  and  is  consequently  called  the  “implicit”  portion  of  the 
equation. 

The  spatial  derivatives  in  (20)  are  to  be  approximated  by  appropriate  finite-difference  quotients.  When  this  is  done, 
the  symbol  Un(r)  is  replaced  by  U"  where 

i -i 

r,=  2  Ar,  +  r0 

/  =  i 

(r0  is  the  base  radius,  or  the  radius  of  the  Sun),  A r,  =  ri+ ,  —  r(,  and  r,  =  r0.  In  constructing  these  differences,  it  will  be 
important  to  allow  for  variable  grid  size  (a  slowly  varying  time  step  is  already  permitted  with  the  scheme  invoked, 
while  retaining  second  order  accuracy).  For  the  computations  used  here,  the  following  three-point  second-order 
accurate  central-difference  approximations  are  used: 

(If), 4  xr+/'(s^T_^)_fc]+°[4r'!i'  ,2la) 

Equation  (21b)  results  from  combining  backward  first-derivative  approximations  at  A  r,  and  Ar,+  1  with  forward 
first-derivative  approximations  at  A  r,_!  and  A  r,  to  find  the  second  derivative.  This  gives  a  form  in  which 

can  be  factored  out.  A  more  conventional  derivation  is  to  use  the  average  of  A  r,  and  A  r,  _ , ,  together  with  two-point 
centered  first-derivative  approximations  over  Ar,  and  Ar,_ ,,  to  find  the  second  derivative  on  a  variable  grid  (Schnack 
1978).  These  two  approximations  are  equivalent  to  within  0(\rf)  times  the  rate  of  change  of  Ar,  with  radius.  Both 
(21a)  and  (21b)  are  strictly  second  order  accurate  only  for  slowly  varying  grid  spacing.  At  the  boundaries,  (21a)  and 
(21b)  will  be  impossible  to  use,  so  at  those  grid  points  the  following  one-sided  three-point  backward  and  forward 
(respectively)  approximations  are  used: 


(«)- 


~(2Ar,_,Ar,.2  +  Ar*_2)/  +(Ar,_,  +  Ar,,;)2/,,  -  Ar,2_,/;_2 
[A»)-lArj_2(ArJ_ ,  +  Ar4_2)] 


9£\  _  ~ (2Ar4Ar,-+ ,  +  Ar^ ,)/,+( Ar,+  Ar^,)2/^,  —  Ar,2/,42 
drj,  [  Ar,Ar1  + ,(  Ar,  +  Arl4. ,)] 


r 
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where  the  backward  (forward)  differencing  scheme  is  to  be  applied  at  the  outer  (inner)  grid  point.  Second  derivatives  at 
the  boundaries,  required  when  including  dissipation,  will  be  discussed  in  §  IV. 

Use  of  (21)  in  (20)  produces  a  block-tridiagonal  system  of  equations  of  the  form: 

xr  ,Af ,  +  Y”MJ"  +  ZA  ,A V,\  ,  =  T,\  (23) 

where  X"_ , ,  Y,",  and  Z,"  ,,  are  matrices  of  dimension  3X3,  and  the  T"  are  vectors.  This  system  of  equations  can  be 
solved  using  standard  algorithms  such  as  that  described  by  Isaacson  and  Keller  (1966)  or  by  Varga  (1962).  In  this 
application,  the  BW  formulation  contains  three  levels  of  data  ( n  +  1, «,  and  n  —  1);  but  only  two  levels  of  data,  1)"  and 
At/”,  need  be  stored  for  each  spatial  grid  point.  This  decreases  the  amount  of  computer  memory  required  for  these 
matrices  by  one- third. 

At  the  start  of  the  calculation,  it  is  assumed  that  a  complete  description  of  the  solar  wind  flow  is  provided  at  each 
grid  point.  This  provides  the  initial  state,  but  at  only  one  time  level.  Consequently,  as  described  earlier,  X  =  |  and 
i  -  0.0  for  the  first  time  step.  Thereafter,  they  are  taken  to  be  1.0  and  0.5,  respectively.  The  initial  state  is  essentially 
arbitrary — it  need  not  be  a  solution  to  the  steady  solar  wind  equations,  although  it  will  be  such  when  modeling  the 
transients  described  in  §  VI.  Various  choices  will  be  discussed  in  §  IV. 

IV.  ADDING  DISSIPATION 

The  implicit  generalized  time-differencing  scheme  (16)  is,  as  applied  here,  temporally  dissipative  except  for  the 
longest  and  shortest  wavelengths  (BW,  §  V).  Since  the  phase  error  of  the  short  wavelengths  is  large,  it  is  sometimes 
necessary  to  add  dissipative  terms  to  the  explicit  portion  of  the  equations  to  damp  the  short  wavelengths.  Apparently 
this  introduces  a  new  stability  bound  (Steger  1978).  so  Steger  (private  communication)  suggests  that  a  dissipative  term 
also  be  added  to  the  implicit  portion  of  the  equations.  The  terms  chosen  here,  which  are  of  the  same  form  as  those  used 
by  Steger  (1978)  and  Steger  and  Bailey  (1980),  are  shown  below  in  the  modified  form  of  equation  (20)  to  be  used 
henceforth  for  the  analysis: 

{/  +  y£|  |;(^+C-/*r)"  +  ^(/*)',+(£)'’-f,Ar,2^(vA),  JaI/" 

=  -~^|^;(F+F)"+(W)"-f«rAr,4^;(vA)2}t/'’  + T^Af/'’*1.  (24) 

The  operators  A  and  v  have  the  conventional  definitions: 

(A)'/=^T’  (25a) 

(25b) 

so  that  they  indicate  forward  and  backward  differences,  respectively.  The  scaling  parameters  r,  and  uc  are  the  radius  of 
the  adiabatic  (or  polytropic)  critical  point  and  the  adiabatic  sound  speed  at  the  critical  point  computed  using  the  base 
temperature  ( T0 )  and  defining  the  critical  point  in  the  manner  used  by  Parker  (1963): 


GM 

4  yRT0 ' 

(26a) 

=  (2y  RT0)'/2. 

(26b) 

These  definitions  are  used  so  that  it  is  not  necessary  to  incorporate  any  additional  scaling  in  the  use  of  e,  and  te. 

Fourth  order  dissipation  terms  have  been  added  to  the  explicit  portion  of  (24).  This  does  not  impair  the  formal 
accuracy  of  the  method.  In  contrast,  second  differences  operating  on  A U"  have  been  added  to  the  implicit  portion  of 
(24).  Implicit  use  of  fourth-order  differences  as  dissipation  terms  would  probably  produce  even  wider  stability  bounds 
(Steger  1978),  but  would  require  block-pentadiagonal  matrix  inversions  of  significantly  greater  complexity.  An  analysis 
of  the  stability  bounds  on  et  and  e,  has  not  been  made  for  this  specific  study.  However,  Steger  (1978)  has  found  that 
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typically  e,  =  0[1]  and  e(  =  2.5  ef  are  upper  bounds,  but  that  both  terms  may  be  much  larger  during  an  impulsive  initial 
start.  For  the  computations  done  here  with  the  solar  wind  equations,  an  attempt  has  been  made  to  include  as  little 
dissipation  as  possible.  Consequently,  values  of  ee  and  e,  have  been  chosen  to  be  less  than  or  equal  to  2.0  X  I  (T 2  for  all 
of  the  results.  Abbreviated  studies  have  shown  no  effects  larger  that  1  %  in  timing  or  magnitude  of  the  profiles  for 
values  of  either  t,  or  e,  up  to  at  least  0(1]. 


V.  BOUNDARY  CONDITIONS 

A  choice  must  be  made  for  the  boundary  conditions  to  be  applied  both  at  the  base  of  the  corona— at  the  first  grid 
point — and  at  the  outer  grid  point  of  the  computational  mesh.  For  all  the  examples  considered  here,  the  radius  of  the 
outer  grid  point  has  been  chosen  such  that  the  flow  there  is  supersonic,  and  the  additional  condition  is  made  that  the 
flow  be  subsonic  at  the  base.  Neither  of  these  conditions  is  necessarily  always  true  for  any  type  of  coronal  transient, 
but  will  be  required  to  be  the  case  here  in  order  to  avoid  switching  the  manner  in  which  the  boundaries  are  treated. 

An  analysis  illustrating  the  mathematically  proper  method  for  specifying  boundary  conditions  for  time-dependent 
solutions  to  the  solar  wind  equations  describing  adiabatic  (or  polytropic)  flow  was  presented  by  Nakagawa  and 
Steinolfson  (1976).  The  application  of  this  and  other  methods  was  described  by  Steinolfson  and  Nakagawa  (1976).  The 
essence  of  the  analysis  is  (see  also  Courant  and  Hilbert  1962)  that  for  a  mixed  initial  boundary  value  problem  only  as 
many  boundary  conditions  can  be  specified  as  there  are  characteristics  issuing  into  the  region  of  interest  from  that 
boundary.  For  the  present  problem,  there  are  three  real  characteristics  which  are  related  to  (i)  the  flow  speed,  (ii)  the 
flow  speed  minus  the  sound  speed,  and  (iii)  the  flow  speed  plus  the  sound  speed.  Since  the  flow  is  subsonic  at  the  base, 
the  flow  speed  minus  the  sound  speed  is  negative,  resulting  in  only  two  rather  than  three  characteristics  issuing  into  the 
region  of  interest.  Thus  only  two  dependent  variables  can  be  specified  at  the  base  of  the  corona.  According  to  the 
theory  of  characteristics,  the  third  dependent  variable  at  the  boundary  must  be  determined  using  some  sort  of 
extrapolation  depending  on  the  flow  variables  interior  to  the  boundary. 

In  the  present  case,  the  dependent  variables  at  the  inner  boundary  can  be  taken  as  the  temperature,  density,  and 
velocity,  and  the  choice  is  made  to  use  temperature  and  density  as  the  specified  boundary  variables.  It  follows  now  that 
the  velocity  at  the  base  must  be  found  using  an  extrapolation  from  the  interior  grid  points. 

Several  methods,  in  addition  to  the  mathematically  “proper  method”  using  the  characteristic  equations,  are 
examined  by  Steinolfson  and  Nakagawa  (1976).  All  of  these  methods  have  been  tested  in  the  present  problem,  although 
the  characteristic  method  was  used  only  for  polytropic  flow.  The  most  reliable  extrapolation  from  the  viewpoint  of 
being  able  to  find  the  smoothest  numerical  solutions  was  simply  to  let 

(mr)"=(mr)"v  (27a) 

where  i  =  1  and  i  =  2  indicate  the  grid  points  at  the  inner  boundary  and  adjacent  to  that  boundary,  respectively,  on  a 
grid  running  from  i  =  1  to  /  =  /.  This  extrapolation  is  exactly  true  for  steady  flow,  but  is  only  a  crude  approximation 
for  time-dependent  flow  that  becomes  more  accurate  as  the  size  of  the  first  grid  step  decreases.  The  reason  this 
extrapolation  works  well  is  that  it  is  the  least  sensitive  to  the  detailed  nature  of  the  profiles  and  involves  only  two  grid 
points.  Higher  order  or  characteristic  method  schemes  require  more  complex  extrapolations,  to  which  the  numerical 
stability  of  implicit  schemes  can  be  sensitive. 

At  the  outer  boundary,  the  flow  is  supersonic.  In  principle,  again  from  characteristic  theory,  no  information  can 
propagate  upstream  because  all  characteristic  speeds  are  positive  while  the  region  of  interest  is  in  the  negative  direction 
from  the  outer  boundary.  Thus,  an  extrapolation  from  the  interior  flow  field  profiles  must  be  made  for  all  three  of  the 
variables.  Here,  the  very  simplest  extrapolation  is  used  (see,  e.g.,  BW)  by  taking 

Up  =  £//L , ,  (27b) 

where  i  =  /  and  i  =  /  - 1  are  the  outer  boundary  and  adjacent  grid  points,  respectively.  This  extrapolation  has  been 
found  here  to  be  extremely  reliable,  but  the  solution  can  depend  on  how  the  temperature  gradient  at  the  outer 
boundary  is  chosen  or  computed  in  finite  difference  form.  This  problem  will  be  discussed  more  in  §  VI. 

Conditions  (27a,b)  can  also  be  used  to  advance  the  boundary  conditions  to  the  next  time  step,  giving 

Wp  =  AUf_it  (28a) 


(im,)i=(^f)2. 
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Because  T0  and  p0,  the  temperature  and  density  at  the  base,  are  specified,  then 
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Afc;-  =  0,  (28c) 

(28d) 

This  completes  the  formal  statement  of  the  boundary  conditions. 

The  application  of  (24),  together  with  (27)  and  (28),  at  the  boundaries  of  the  computation  region  is  best  thought  of  in 
terms  of  stepping  through  the  sequence  of  operations  required  to  advance  the  solution  from  time  level  n  to  level  n  +  1 . 
Consider  first  the  explicit  portion  of  the  boundary  condition,  and  then  the  implicit  portion. 


a )  Explicit  Portion 

The  i  =  1  inner  boundary  is  first  encountered  for  /  =  2  when  inserting  the  spatial  differencing  (2 1  a)  into  the  right 
hand  side  of  (24).  Evaluating  the  derivative  8 (F+  Y)/d r  at  /  =  2  requires  values  for  F"  and  V".  Using  (A2)  and  (13),  it 
is  seen  that  F"  can  be  determined  algebraically  in  combination  with  the  boundary  conditions  (27a,  b).  However,  V" 
requires  an  estimate  of  (8T/8r)7  which,  using  (A6g),  in  turn  requires  an  estimate  of  (8f>/8r),,  (8m,/8r)| ,  and 
(8e/8 r)".  To  find  these,  the  one-sided  differences  in  (22b)  are  used  and  the  required  derivates  at  the  boundary  point 
are  found  using  (b,  mr,  e),  at  /  =  1,  2,  and  3. 

Information  at  the  boundary  is  also  required  for  the  dissipation  term  (vA )]ll".  To  evaluate  this  fourth-order 
difference,  a  centered  five-point  scheme  is  invoked  in  the  application  of  (25a,  b).  At  /  =  2,  this  requires  not  only 
information  at  the  i  =  1  grid  point,  but  also  at  an  imaginary  grid  point  /  =  0.  The  i  ~  1  information  is  found  in  the  same 
way  as  above.  The  approach  used  here  for  the  /  =  0  point  is  to  set  U0n  —  U"  or,  equivalently,  to  take  (v)U"  —  0. 
Similarly,  at  /  =  /  —  1  the  centered  five-point  difference  used  in  evaluating  the  dissipation  term  requires  information  at 
the  imaginary  grid  point  i  =  /  +  1 .  Precisely  the  same  extrapolation  is  used  here  as  for  i  =  2  such  that  U,\  ,  =  U,"  or 
(A  )V,"+]  =  0. 


b )  Implicit  Portion 

In  application  of  (21  a, b)  to  the  left  hand  side  of  (24),  the  three-point  central  difference  scheme  results  in  a  now 
spatially  differenced  equation  of  the  form  (see  also  [23]): 


,At/;_ ,  +  y; A(/;  +  Z/V  .a ur+  ,  =  r,\ 
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(33) 

where  /=  the  3X3  identity  matrix,  and  e,  is  the  implicit  dissipation  coefficient  in  (24)  and  is  unrelated  to  the  /-index  of 
the  grid  point. 
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For  i  =  2  and  i  =  I  - 1,  the  boundaries  are  encountered  in  the  terms  X". ,  A U"  ,  and  Z/V  ,Al/("  „  respectively.  At  the 
outer  boundary,  A U?  -  A 1//L ,  from  (28a).  At  the  inner  boundary,  the  treatment  is  slightly  more  complex.  A U"  is  found 
from  (28b),  (28c),  and  (28d): 


At/"  = 


To  incorporate  this  into  (29)  it  is  necessary  to  write  X"AU”  in  the  form  X"K"„AU".  That  is: 

K"BAU2  =  At/,". 


This  is  achieved  if 


such  that 


K" 

n  IB 
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0  1  0 
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(34) 


(35) 


(36) 


(37) 


Substitution  of  (34)-(37)  into  (29)  produces  a  system  of  linear  equations  for  A U",  2  <  /  <  /  -  1 ,  which,  in  matrix  form, 
is 


'(Y2  +  XtKlB)  Z3  0 
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'  A L\ 
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(38a) 


or 

«"At/"  =  r\ 


(38b) 


where  At/"  and  T "  are  vectors  of  dimension  ( /  -  2)  and  ®"  is  a  block  tridiagonal  matrix  with  blocks  of  dimension  3X3. 

The  procedure  is  to  compute  fi"  and  T”  for  the  present  time  step  and  to  then  solve  (38)  for  At/".  This  gives  the 
increment  in  the  depe.ident  variable  for  that  time  step.  The  new  value  of  the  dependent  variable  at  the  advance  time 
step  is  then  computed  from  t/"+ 1  -U”  +  At/",  placing  the  old  value  l/"  in  an  appropriate  mass-storage  device  and 
replacing  it  with  t/"+ 1  so  as  to  require  core  memory  only  for  At/"  and  t/"+ '.  Algorithms  for  the  solution  of  (38)  have 
been  described  in  many  texts  such  as  those  by  Isaacson  and  Keller  (1966,  pp.  58-61)  and  Varga  (1962,  pp.  99-102). 


VI.  NUMERICAL  RESULTS  COMPARING  POLYTROPIC  AND  THERMALLY  CONDUCTIVE  FLOW 

The  main  example  will  be  a  comparison  of  a  transient  for  thermally  conductive  flow  with  a  similar  transient  for 
polytropic  flow.  However,  before  this  can  be  done,  the  appropriate  steady  states  (that  will  then  be  used  as  initial  states 
for  the  transient)  must  be  found.  This  first  step  is  also  done  using  the  implicit  algorithm  in  order  to  demonstrate  its  use 
in  finding  steady-state  solutions  to  the  solar  wind  equations. 
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In  these  examples,  the  mass  of  the  Sun  is  1.989X  1033  g,  and  the  inner  boundary  radius  is  6.96X  10'°  cm.  Also,  the 
spatial  grid  was  common  for  all  cases,  being  described  by  the  formula: 


A  r, .  ,  =  Ar,  1  + 


^pc*  \ 

ioo  r 


7-1>i>1, 


(39) 


for  A  r,  <  A  rmax ,  and  A  r,  —  A  rmax  at  larger  radii.  In  (39),  A  r,  is  specified,  as  is  A  rpct ,  the  percent  that  A  r,  is  allowed  to 
increase  at  each  grid  step  up  to  its  also  specified  maximum  value  A  rmax.  The  maximum  radius,  rmax,  is  prescribed 
before  the  grid  is  computed.  Normally,  the  actual  grid  would  exceed  rmax  by  some  fraction  of  A  rmax  if  so  allowed.  Here, 
the  grid  is  terminated  at  the  largest  radius  that  is  less  than  rmax  by  the  latest  grid  step  size.  For  all  the  examples, 
Ar,  =  0.01  Rq,  Arpc[  =  10%,  rmM  =  26  R0,  and  A =  2.0  Ra.  This  results  in  a  grid  of  59  points  spanning  1.00  R„  to 
25.697  Rq,  with  the  last  three  steps  being  the  maximum  size  of  2.0  Re.  This  gives  good  resolution  near  the  base  of  the 
grid,  with  only  1 1  grid  points  lying  between  10  RQ  and  26  RC). 

In  the  first  examples  of  relaxation  from  an  initial  state  to  the  steady  state,  the  initial  state  will  be  chosen  as 
isothermal  flow  for  both  the  polytrope  and  with  thermal  conduction.  The  base  temperature  and  density  are  1.6X  10fc  K 
and  7.4 X  107  cm-3,  respectively— the  same  values  chosen  by  Whang  and  Chang  (1965)  in  one  of  the  earliest  solutions 
to  the  steady  state  solar  wind  equations  with  thermal  conduction.  As  described  in  §  II,  the  thermal  conduction 
description  used  is  the  collision-dominated  one  given  by  Spitzer  (1962)  for  a  fully  ionized  hydrogen  plasma.  The 
poly  tropic  index,  gamma,  is  taken  to  be  1.10  specifically  because  this  will  produce  a  steady  state  very  nearly  the  same 
as  that  with  thermal  conduction,  between  1.0  and  21.0  RQ.  (If  a  much  larger  index  were  used  or  required,  a  Parker-type 
solution  would  not  exist;  see,  e.g.,  Parker  1963,  pp.  59ff.) 

In  these  relaxations,  the  initial  state  is  essentially  arbitrary,  under  the  restrictions  stated  in  §  V.  The  boundary 
conditions  are  held  constant,  but  the  physics  is  changed  at  the  initial  instant.  Hence,  little  physical  significance  can  be 
attached  to  the  intermediate  steps  between  the  initial  state  and  the  steady  state.  The  only  significance  that  exists  is  that 
the  amount  of  physical  time  taken  to  reach  the  steady  state  is  a  reflection  of  the  true  amount  of  time  the  corona  needs, 
under  the  prescribed  physics,  to  relax  to  a  steady-state  after  some  perturbation.  To  minimize  the  computational  time 
needed  to  reach  the  steady  state,  the  time  step  is  allowed  to  expand  in  a  manner  similar  to  the  grid  size  in  (39),  except 
that  the  time  step  is  allowed  to  grow  only  if  the  profiles  are  changing  no  more  rapidly  than  some  specified  amount 
anywhere  on  the  grid.  In  general,  the  steady  states  are  achieved  in  physical  times  of  less  than  30  hours  out  to  20  Ra. 
This  has  been  found  to  require  about  18  times  steps  of  up  to  2.0  hours  length  in  physical  time.  Comparing  this  with 
solutions  of  the  steady  state  equations,  it  is  analogous  to  the  number  of  iterations  necessary  on  the  location  of  the 
critical  point  in  order  to  match  a  given  temperature  and  density  at  the  base  of  the  corona.  The  computational  effort  for 
this  relatively  simple  set  of  equations  is  very  nearly  the  same  in  either  case  (Kopp  and  Steinolfson,  private 
communication). 

The  results  of  the  two  relaxations  are  shown  in  Figure  1,  giving  the  profiles  out  to  21  R0.  This  is  done  to  absolutely 
avoid  any  contamination  from  the  treatment  of  the  outer  boundary  condition.  Curve  (i),  in  each  panel,  gives  the  initial 
state  temperature,  velocity,  and  density,  respectively.  Curves  (ii)  and  (iii)  show  the  resultant  steady  states  for  thermally 
conductive  flow  and  polytropic  flow  (gamma  =  1.10),  respectively,  after  a  physical  time  of  40  hours.  As  must  be  the 
case,  these  profiles  match  the  known  solutions  for  these  boundary  conditions,  except  for  small  differences  in  curves  (ii) 
due  to  the  retention  of  the  Coulomb  logarithm  in  the  description  of  thermal  conduction,  when  comparing  with  the 
results  of  Whang  and  Chang  (1965). 

It  is  of  some  interest  from  a  numerical  point  of  view  that  these  results  are  produced  with  no  knowledge  of  the 
asymptotic  nature  of  the  flow  profiles  at  large  distances  from  the  Sun  and  with  no  reference  to  critical  point  nature.  It 
is  also  notable  that  no  effort  was  made  to  choose  an  initial  state  near  the  final  state  in  the  relaxations.  The  steady  states 
were  found  easily  from  a  distant  initial  state  with  time  steps  that  correspond  to  Courant  numbers  of  up  to  0[1O2-1O3] 
even  for  polytropic  flow.  As  a  test  of  their  stability,  it  is  also  possible  to  begin  with  either  curve  (ii)  or  (iii)  as  an  initial 
state,  appropriately  change  the  physics,  and  relax  to  curve  (iii)  or  (ii),  respectively. 

Curves  (ii)  and  (iii)  can  now  be  used  as  initial  states  for  the  study  of  transients.  The  example  that  will  be  shown  here 
invokes  an  increase  in  the  base  temperature  from  1.6  x  106  K  to  2.0  X  lO6  K  at  /  =  0.  and  then  holds  the  density  and 
temperature  constant  thereafter.  This  impulsive  change  in  the  base  temperature  is  not  large  enough  to  produce  a  shock, 
so  a  smooth  transient  in  the  flow  will  propagate  outwards.  As  opposed  to  the  above  relaxation  to  find  the  steady  states, 
the  physical  description  is  now  constant  but  the  boundary  conditions  are  changed— producing  a  transient  whose 
properties  are  of  physical  interest.  The  results  of  this  exercise  are  shown  in  Figure  2.  Figure  2  a  shows  the  transient 
induced  on  polytropic  flow,  and  Figure  2  b  shows  the  same  for  thermally  conductive  flow.  Profiles  are  plotted  at  t  =  0, 
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Fig.  I.—  (i)  Isothermal  solar  wind  flow,  used  as  the  initial  state  for  the  relaxations  to  steady  states  with  constant  boundary  conditions 
shown  as  curves  (til  and  (iii).  (it)  Final  state  (or  thermally  conductive  flow,  (iii)  Final  stale  for  polvtropic  flow  (gamma  =1.10).  Both  final 
states  are  for  relaxations  times  of  slightly  more  than  40  hours. 


2,  8,  and  20  hours  in  each  case.  It  is  immediately  obvious  that  the  two  examples  produce  quite  different  results  even 
though  they  begin  with  similar  initial  states.  Most  of  these  differences  are  due  to  the  ability  of  a  thermal  pulse  to 
propagate  at  many  times  the  sound  speed.  However,  some  of  the  differences  are  due  to  the  specific  values  chosen  for 
the  base  temperature  and  density,  as  will  be  explained  in  more  detail  shortly. 

The  transient  in  Figure  2 a  is,  except  for  its  small  amplitude,  similar  to  those  in,  e.g.,  the  model  described  by 
Steinolfson  and  Nakagawa  (1976)  The  initial  rise  in  temperature  produces  a  velocity  pulse  and  compression  which  is 
evident  in  the  density  profiles.  Tin.  initial  velocity,  temperature,  and  density  increases  are  all  in  phase  because  no  signal 
propagates  more  rapidly  than  the  sound  speed.  This  transient  is  also  very  much  like  transients  produced  in  models 
limited  only  to  regions  of  supersonic,  polytropic  flow.  In  other  words,  there  is  nothing  new  here,  except  in  comparison 
between  Figures  2 a  and  2b. 

In  Figure  26,  it  is  obvious  that  the  high  thermal  conductivity  allows  a  thermal  pulse  to  run  out  in  front  of  the 
velocity  and  density  transients.  In  fact,  this  effect  is  highly  exaggerated  because  the  collisional  thermal  conductivity 
used  here  is  not  a  realistic  description  of  the  mid-  and  outer  corona.  However,  this  is  not  considered  to  be  a  deficiency 
because  the  primary  intent  herein  is  to  demonstrate  the  technique  of  including  thermal  conduction  and  to  illustrate 
comparison  of  this  classical  case  with  polytropic  flow.  In  .fie  transient  with  thermal  conduction,  there  is  an  initial 
temperature  rise  at  all  radii  after  just  a  few  minutes.  This  temperature  rise  is  sr  n  ihe  plots  of  temperature  versus 
time  at  several  different  radii  in  Figure  3.  Subsequently,  the  increased  tempera  "tuses  a  velocity  transient  to  be 
initiated  at  the  base,  much  as  with  the  polytrope.  This  is  shown  in  the  velocity  pronie  at  r  time  of  2  hours.  At  the  front 
of  the  velocity  transient  is  a  compression  region  which  locally  increases  the  temperature  Jowever,  unlike  in  polytropic 
flow,  this  local  temperature  increase  is  quickly  communicated  outward  to  essentially  i.u  '-.ii.  Behind  the  compression 
is  a  rarefaction  that  produces  a  local  temperature  minimum.  In  front  of  the  compression  the  now  enhanced 
temperature  causes  plasma  to  begin  moving  at  all  radii— a  phenomenon  which  cannot  occur  for  a  polytrope,  and  which 
corresponds  to  conductive  damping  of  the  principal  velocity  pulse.  Thermal  energy  has  been  conducted  outward  from 
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Fig.  2 —  The  passage  of  a  transient  on  (a)  polvtropic  flow  (gamma  =  1.10)  and  (/>)  thermally  conductive  flow.  The  initial  states  arc  those 
shown  as  curves  (ii)  and  (iii)  in  Fig.  I.  Profiles  arc  shown  at  elapsed  time  of  0.  2.  8.  and  20  hours  in  all  cases  The  transient  was  introduced  by 
instantaneously  increasing  the  base  temperature  from  1.6  x  10*  to  2.0  x  10*  K  t  =  0.  and  holding  the  temperature  and  density  at  the  base 
constant  at  all  later  times.  Unlike  Fig.  2a.  the  final  state  in  Fig.  2b  has  a  flow  velocity  lower  than  the  initial  state.  This  effect  is  due  to  the 
particular  choice  of  base  temperatures  and  density  for  the  thermally  conductive  flow,  and  enhances  the  apparent  differences  between  Figs.  2a 
and  2  b. 
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Fig.  3.  — The  time  variation  of  the  temperature  at  1.5.  2.5,  5.0,  10.0.  and  20.0  solar  radii  for  the  example  of  a  transient  in  thermally 
conductive  flow  shown  in  Fig.  2b. 
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the  principal  velocity  pulse  and  is  being  partially  converted  to  bulk  motion.  The  density  transient  is  similar  to  the 
transient  for  polytropic  flow,  although  the  differences  are  deemphasized  in  this  logarithmic  plot. 

At  later  times  the  differences  between  Figure  3a  and  3 b  are  just  an  extension  of  the  above  description.  In  addition, 
the  pathological  nature  of  the  initial  choice  of  boundary  conditions  for  thermally  conductive  flow  becomes  evident.  In 
Figure  3 b,  the  velocity  at  20  hours  is  seen  to  be  lower  than  the  initial  velocity.  This  is  a  real  effect  in  that  for  the 
specified  base  temperature  and  density,  a  slight  increase  in  the  base  temperature  will  result  in  a  lower  velocity  at  1  AU 
and  at  all  radii  out  to  1  AU,  as  described  by  Holzer  and  Leer  (1981).  This  is  radically  different  from  the  case  for 
polytropic  flow,  which  produces  increased  velocities  for  increased  base  temperatures  in  all  cases.  A  significant  portion 
of  the  difference  in  velocity  profiles  between  Figures  3a  and  3 b  must  be  attributed  to  this  effect. 

This  pathological  effect  can  be  suppressed  by  choosing  initial  and  final  base  temperatures  and  densities  that  result  in 
similar  initial  and  final  profiles  for  both  polytropic  and  thermally  conductive  flow.  Because  the  problem  is  highly 
nonlinear,  this  suppression  is  not  entirely  possible.  However,  a  good  approximation  to  such  a  case  is  shown  in  Figure  4. 
Here,  the  base  density  has  been  chosen  to  be  a  lower  value  of  1.0X107  cm  \  while  the  initial  and  final  base 
temperatures  are  the  same  as  in  Figure  2.  Inspection  of  Figure  1  from  Holzer  and  Leer  (1981)  shows  this  choice  of 
parameters  to  result  in  increasing  velocities  with  increasing  base  temperature  and  thermal  conduction.  The  initial  and 
final  profiles  are  quite  similar,  although  not  identical.  The  results  for  a  polytrope  are  shown  in  Figure  4a  and  for 
thermally  conductive  flow  in  Figure  4b.  Profiles  are  plotted  at  t  =  0,  2,  5,  and  10  hours  in  Figure  4b.  In  Figure  4a,  an 
additional  profile  at  t  =  20  hours  is  plotted  because  it  was  found  that  it  took  longer  for  the  polytrope  flow  to  reach  an 
apparent  steady  state.  Although  the  time  of  arrival  of  the  transient  is  slightly  different,  the  qualitative  nature  of  the 
transient  on  polytropic  flow  is  the  same  as  in  Figure  2.  For  conductive  flow,  there  are  some  large  differences  between 
this  case  and  the  case  shown  in  Figure  2 — as  was  expected.  The  two  most  striking  differences  are  that  the  magnitude  of 


Fig.  4. —  The  passage  of  a  transient  on  (a)  polytropic  flow  (gamma  =  1.05)  and  ( b)  thcrmalK  conductive  flow.  The  initial  states  have  the 
same  base  temperature  as  curve  (ii)  and  (iii)  in  Fig.  I.  but  the  base  densities  are  now  1.0*  It)1  cm  '  rather  than  7  4  *  It)7  cm  '.  Profiles  arc 
shown  at  elapsed  times  of  0,  2,  5.  and  10  hours  in  Fig.  4 b,  and  at  the  same  times  plus  at  20  hours  in  Fig  4a  because  the  polytropic  flow 
relaxes  more  slowly  to  a  steady  state  than  does  the  thermally  conductive  flow  The  transient  was  again  introduced  by  instantaneously 
increasing  the  base  temperature  from  1.6X  I06  to  2.0X  10'’  K  atr  =  0.  and  holding  the  temperature  and  density  at  the  base  constant  at  all  late 
times.  In  addition  to  the  initial  states,  the  final  state  in  Fig.  4u  is  now  similar  to  the  final  state  in  Fig.  4b  The  differences  between  the  two 
examples  are  primarily  due  to  the  ability  of  a  thermal  pulse  to  propagate  out  in  front  of  the  velocity  pulse  in  Fig.  4 b 
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the  temperature  rise  in  front  of  the  compression  region  is  much  larger  and  that  the  velocity  pulse  is  much  smoother  and 
more  spread  out.  These  differences  are  due  to  the  overall  lower  density  in  the  corona  and  the  corresponding  higher 
thermal  conductivity,  causing  a  significant  shift  in  the  energy  distribution  and  propagation.  What  is  more  important  is 
that  this  figure  demonstrates  the  important  effects  of  thermal  conduction  on  the  propagation  of  a  transient  in  the 
corona,  in  comparison  to  a  similar  transient  in  a  polytropic  gas.  The  initial  and  final  steady  states  are  very  close  for  the 
two  opposing  cases,  so  the  differences  represent  a  realistic  picture  of  these  effects. 

It  is  possible  to  view  the  results  shown  in  Figure  4 b  in  a  slightly  different  manner.  Because  the  heat  pulse  propagates 
outward  in  front  of  the  velocity  transient  at  the  expense  of  the  energy  in  the  velocity  transient,  the  net  effect  is 
extraction  of  energy  from  some  instantaneous  state  low  in  the  corona  and  deposition  of  that  energy  in  the  mid-  and 
outer  corona.  In  a  sense,  the  passage  of  the  velocity  pulse  has  a  refrigerating  effect  on  the  low  corona.  This  is  only  a 
transient  effect,  so  it  does  not  necessarily  represent  any  average  heating  of  the  outer  corona  or  cooling  of  the  inner 
corona,  although  such  could  be  the  case  in  the  presence  of  a  continuous  wave  train. 

The  last  two  examples  that  will  be  shown  are  the  result  of  an  investigation  of  the  effect  of  alternate  treatment  of  the 
outer  boundary  condition.  The  results  shown  in  Figures  1-4  used  the  condition  given  in  (27b).  From  (14),  this  would 
seem  to  suggest  that  (3r/9r);,  the  temperature  gradient  at  the  last  grid  point,  is  zero.  However,  this  is  not  necessarily 
the  case,  depending  on  how  derivatives  are  constructed  in  finite  difference  form.  In  the  previous  examples,  derivatives 
at  /  =  /  were  constructed  using  the  three-point  one-sided  form  given  in  (22a)  which,  when  combined  with  (27b),  results 
in  a  value  for  the  temperature  gradient  at  the  outer  grid  point  that  is  approximately  half  the  gradient  at  the  next  point 
inward.  That  is. 


The  general  requirement,  physically,  is  that  the  outer  boundary  condition  on  the  temperature  gradient  should  allow 


Flo.  5  — At  attempt  to  find  a  steady  state  from  an  initial  state  of  identical  to  curve  (i)  in  Fig.  I.  Here,  however,  the  boundary  condition  at 
the  outer  boundary  was  altered  to  make  the  temperature  gradient  as  computed  between  the  /  =  /  and  i  =  (/  —  !)  grid  points  equal  that 
between  the  i  =  ( /  -  1)  and  i  =  (  /  -  2)  grid  points.  Curves  (i)  show  the  initial  stat",  and  curves  (ii)  show  the  profiles  after  23  hours.  These 
profiles  continue  to  diverge  from  the  known  final  steady  state  with  time. 
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Fig.  6. — Temperature  variation  for  a  transient  on  thermally  conductive  flow  two  hours  after  the  temperature  at  the  base  was  impulsively 
increased  from  1.6X  106  to  2.0x  I06  K.  The  base  density  is  7.4X  I07  cm  \  Two  initial  states,  curves  (i)  and  (ii).  have  been  used.  Curve  (i)  is 
the  same  initial  state  shown  in  Fig.  1  b.  Curve  (ii)  is  the  same  as  (i)  except  with  outer  boundary  condition  as  in  Fig.  5.  Here,  a  relaxation,  using 
curve  (i)  as  an  initial  state,  did  produce  the  new  steady  state  (ii)  with  slightly  lower  temperatures  near  the  outer  boundary.  The  two  transients 
seem  similar,  in  spile  of  the  modified  boundary  condition.  This  will  not  be  the  case  as  the  transient  approaches  the  boundary  of  the 
computation  region. 

heat  to  conduct  through  the  boundary  in  as  natural  a  manner  as  possible.  Alternate  treatment  of  the  finite  differencing 
at  /  =  /  can  greatly  affect  how  well  this  requirement  is  met,  regardless  of  the  physical  treatment  of  the  boundary 
condition.  It  can  be  shown  analytically  that  nearby  values  of  37"/ 9r  to  that  given  in  (40)  do  not  greatly  affect  the 
solution.  The  largest  effects  are  found,  in  any  case,  near  the  boundary  and  when  variables  are  rapidly  changing  near  the 
boundary.  What  cannot  be  done  is  to  change  the  sign  of  either  the  slope  of  the  temperature  gradient  as  computed  using 
the  points  i  =  /  and  t  =  /  —  1 ,  or  of  the  curvature  as  computed  using  the  points  /  =  /,/  —  1 ,  and  1  —  2.  Thus,  choosing  to 
evaluate  dT/dr  using  only  i  —  I  and  i  —  I  —  1,  together  with  the  boundary  condition  (27b),  would  usually  cause 
erroneous  results  because  the  slope  has  been  changed  from  a  negative  value  to  zero.  Similarly,  choosing  (dT/dr),= 
(37y  3r)r_ ,  would  cause  the  curvature  at  the  outer  grid  point  to  be  computed  as  zero.  The  results  of  this  second  choice 
are  what  are  shown  in  Figures  5  and  6. 

The  choice  that  (dT/dr),  =  (dT/dr),_^  might  be  expected  to  be  a  reasonable  alternative  to  the  extrapolation  (40). 
In  fact,  an  extrapolation  similar  to  (27b)  could  be  invoked  which  would  give  exactly  this  temperature  gradient  when 
applying  (22a).  Physically,  however,  this  choice  represents  an  extreme  enhancement  of  the  manner  in  which  heat  is 
conducting  through  the  outer  boundary.  If  the  temperature  gradient  is  too  steep,  then  (i)  more  heat  will  conduct  out 
past  the  last  grid  point  than  should,  (ii)  temperatures  in  the  vicinity  of  the  boundary  will  be  too  low,  and  (ii)  velocities 
will  be  locally  increased.  Conversely  (and  by  far  a  worse  situation),  if  the  temperature  gradient  at  the  outer  boundary 
becomes  positive  during  a  transient,  then  the  boundary  appears  as  a  heat  source  under  the  above  extrapolation  of  the 
temperature  gradient.  This  heat  source  at  the  outer  boundary  apparently  then  never  disappears  at  later  times.  Such  a 
case  is  shown  in  Figure  5,  where  a  relaxation  from  isothermal  flow  to  a  steady  state  was  attempted  using  precisely  the 
same  parameters  as  in  Figure  1  for  a  relaxation  with  thermal  conduction  (relaxation  from  profiles  [i]  to  [ii]  in  Fig.  1).  In 
Figure  5,  the  initial  isothermal  flow  is  given  by  curves  (i),  and  the  profiles  after  23  hours  physical  time  are  shown  by 
curves  (ii).  The  temperature  has  risen  above  the  isothermal  state,  and  is  continuing  to  rise  at  all  radii  because  at  the 
initial  time  step  the  flow  was  hotter  than  in  the  steady  state  and  so  was  compressed  slightly  near  the  outer  boundary. 
This  compression  raised  the  temperature  slightly  at  the  outer  boundary,  causing  the  temperature  gradient  there  to 
become  positive.  The  consequence  is  as  described  above.  The  outer  boundary  now  acts  as  a  heat  source,  the 
temperature  and  velocity  are  too  large,  and  a  steady  state  is  not  being  approached. 

A  contrasting  situation  is  shown  in  Figure  6.  Here  only  the  temperature  is  plotted  for  two  different  cases,  at  time 
t  =  0  and  /  -  •  2  hours.  The  initial  state  (i)  is  the  same  as  in  Figure  2b.  This  initial  state  was  then  used  in  two  ways.  First, 
a  transient  was  introduced  by  increasing  the  base  temperature  from  1.6X  106  K  to  2.0X  106  K,  with  the  resulting 
temperature  profile  after  2  hours  being  shown  by  curve  (iii).  Next,  (i)  was  used  to  generate  a  new  steady  state  with  the 
alternate  extrapolation  for  the  temperature  gradient  at  the  outer  boundary  that  was  used  for  the  results  shown  in 
Figure  5.  In  the  present  case,  the  temperature  gradient  at  i  =  I  never  becomes  positive,  and  a  new  steady  state  is 
achieved  with  slightly  depressed  temperatures  near  the  outer  boundary— shown  by  curve  i  Then  curve  (ii),  and  the 
associated  velocity  and  density,  were  used  as  the  initial  state  for  another  transient  that  was  introduced  again  by 
increasing  the  temperature  from  1.6X  106  to  2.0  X  106  K.  The  temperature  profiles  after  2  hours  are  shown  by  curve 
(iv),  which  is  essentially  identical  to  curve  (ii)  out  to  3.5  R0.  This  example  shows  that  the  extreme  choice  for 
extrapolating  the  temperature  gradient  at  the  outer  boundary  does  not  greatly  affect  flow  properties  as  long  as  the 
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temperature  gradient  does  not  become  positive  at  that  boundary.  However,  eventually  the  transient  shown  by  curve  (iv) 
will  reach  the  boundary,  causing  a  reversal  in  the  sign  of  the  temperature  gradient.  Subsequent  evolution  of  the 
transient  with  this  boundary  condition  will  continue  to  diverge  from  the  physical  solution  due  to  the  same  effect 
demonstrated  in  Figure  5. 

These  last  two  examples  were  meant  to  demonstrate  the  physical  limitations  on  extrapolation  such  as  (27b)  in 
combination  with  finite  difference  approximations  at  the  boundaries.  As  long  as  these  physical  limitations  are  not 
violated,  it  is  possible  to  show  analytically  that  other  extrapolations  of  the  temperature  gradient  than  that  chosen  here 
only  have  small  effect  on  the  solution,  and  then  only  in  the  vicinity  of  the  boundary. 


VII.  USE  OF  THIS  ALGORITHM  IN  SOLAR  AND  STELLAR  WIND  STUDIES 

Many  processes  can  already  be  modeled  with  the  algorithm  described  in  §§  III-V,  including  propagation  and 
conductive  damping  of  waves  and  wave  trains  that  are  any  combination  of  temperature,  velocity,  and  density  variation 
at  the  base.  In  addition,  alternate  forms  for  thermal  conduction  in  the  context  of  single  fluid  equations  can  be 
investigated,  including  use  of  a  term  derived  from  a  table  look-up  as  opposed  to  a  purely  analytical  form.  However, 
there  are  important  and  obvious  deficiencies  in  the  use  of  equations  ( 1)— (3)  to  describe  solar  and  stellar  winds.  The 
most  important  stellar  winds  are  radiatively  driven,  and  will  be  discussed  below.  For  the  solar  wind,  several  additions 
to  these  equations  might  be  of  interest,  so  a  brief  outline  of  how  to  go  about  some  such  additions  will  be  given  here. 

Many  of  the  interesting  additions  to  the  solar  wind  equations  take  the  form  of  source  (or  loss)  terms  in  the  energy 
and  momentum  equations.  These  terms  are  a  consequence  of  the  inference  that  energy  and  momentum  addition  are 
important  in  the  corona  above  the  transition  region  (Suess  etal.  1977;  Holzer  1977).  It  is  uncertain  what  dependence 
these  terms  have  on  the  conserved  variables,  so  here  arbitrary  functions  will  be  assumed.  These  additions  would  have 
the  effect  of  modifying  the  //-vector  (eqs.  [15]  and  [A4])  to  give 


0 


H  = 


-2  rp  +  bGM/r 2  -  r2qw 
mrGM/r 2  -  r2qe 


(41) 


So  far,  this  does  not  complicate  things  much,  but  it  does  further  remove  the  system  of  equations  from  purely 
conservation  form.  The  additional  complexity  comes  in  evaluation  of  the  nonlinearities  in  H  for  evaluation  at  the 
advance  time  step,  or  in  construction  of  the  Jacobian  matrix  E  (eq.  [19c]).  In  addition  to  £,  there  must  be  one  new 
Jacobian  matrix,  3C,  because  qw  and  qe  may  depend  on  Ur.  This  introduces  the  following  two  matrices: 


TT-3" 
x  dUr ' 

(42a) 

II 

* 

(42b) 

into  appropriate  places  in  equation  (20)  or  (24). 

The  additional  terms  shown  in  (41)  and  (42)  in  solar  wind  studies  suggest  many  new  avenues  of  research  One 
example  will  be  suggested  involving  a  question  of  the  stability  of  a  solution  to  the  equivalent  steady-state  equations. 
Finding  steady  state  solutions  with  totally  ad  hoc  terms  such  as  qw  and  qe  has  been  done  by  Leer  and  Holzer  (1981), 
but  (as  they  point  out)  certain  solutions  involving  an  inferred  standing  shock  transition  from  subsonic  to  supersonic 
flow  in  the  fluid  frame,  in  the  presence  of  extended  heat  and  momentum  addition,  may  not  be  steady.  This  stability 
question  can  be  answered  through  a  relaxation  analysis  such  as  described  in  connection  with  Figure  1,  with  the  above 
additions  to  the  equations. 

The  other  important  addition  to  the  original  equations  that  might  be  suggested  for  one-dimensional  solar  wind 
studies  is  the  extension  to  two-fluid  or  multiple-fluid  flow.  This  requires  additional  equations  for  the  additional  ion 
species.  The  algebraic  complexity  is  somewhat  greater  for  such  models  and  the  dimensions  of  the  vectors  and  matrices 
all  increase  by  at  least  one  for  each  additional  species. 

For  stellar  winds,  the  driving  force  for  those  that  can  be  directly  observed  is  not  the  pressure  gradient,  as  with  the 
solar  wind,  but  rather  radiation  pressure.  This  force  can  be  parameterized  as  being  approximately  proportional  to  the 


local  velocity  gradient  (Castor,  Abbott,  and  Klein  1975),  and  thus  incorporated  into  a  formulation  like  that  outlined 
above  in  (41)  and  (42).  This  allows  direct  numerical  modeling  of  steady  and  transient  flow  in  radiatively  driven  winds. 
Furthermore,  the  functional  relationships  in  (41)  and  (42)  need  not  be  analytic.  Just  as  with  thermal  conduction,  it  is 
completely  possible  to  use  a  table  look-up  to  determine  the  radiative  force  for  any  combination  of  the  variables.  This, 
in  fact,  is  probably  the  preferred  approach  for  stellar  winds  as  the  parameterization  of  Castor,  Abbott,  and  Klein  is 
valid  for  only  a  relatively  small  region  of  parameter  space  (Abbott,  personal  communication). 

In  principle,  generalization  of  the  boundary  conditions  described  in  §  V  would  allow  an  even  broader  treatment  of 
winds.  By  initially  specifying  a  hydrostatic  stratification  and  dropping  the  pressure  at  the  outer  boundary  to  a  very  low 
value  at  /  =  0,  it  should  be  possible  to  study  stellar  wind  startup.  This  has  been  done  already  for  isothermal  flow 
(Kopp,  personal  communication)  and  is  a  problem  of  great  interest.  The  opposite  problem  of  stellar  wind  quenching, 
although  probably  of  less  interest,  could  be  addressed  in  an  analogous  manner. 

The  intent  in  this  section  has  been  to  outline  the  mechanics  of  applying  the  algorithm  described  in  §§  III— V  to 
different  problems.  Other  examples  than  those  given  above  can  be  imagined,  but  most  of  the  necessary  tools  for  any 
possible  example  are  included  in  the  above  discussion.  This  does  not  include  extension  to  more  than  one  spatial 
dimension.  However,  multidimensional  flow  is  discussed  in  great  detail  by  BW. 

VIII.  SUMMARY 

The  time-dependent  solar  wind  equations,  including  thermal  conduction,  present  a  particularly  difficult  problem  in 
their  numerical  solution  due  their  being  extremely  “stiff’— resulting  in  the  Courant  stability  condition  present  with  all 
explicit  time-difference  schemes  requiring  such  small  time  steps  for  numerical  stability  that  it  is  generally  not  feasible 
to  include  thermal  conduction.  Here  a  detailed  description  has  been  given  of  a  general  implicit  time-difference  scheme 
that  completely  avoids  the  above  stability  condition.  The  algorithm  is  mathematically  concise  and  is  found  to  be 
computationally  efficient.  Other  than  the  algebraic  complexity  present  in  all  implicit  schemes,  the  primary  area  in 
which  problems  might  occur  with  this  algorithm  is  in  the  treatment  of  the  boundary  conditions.  A  special  effort  has 
been  made  to  identify  particular  problems  and  solutions  in  dealing  with  the  boundary  conditions  for  varying  cases. 

A  specific  model  is  described  involving  flow  which  is  subsonic  at  the  inner  boundary  and  supersonic  at  the  outer 
boundary  of  the  computation  region.  These  boundary  condition  choices  are  made  because  this  is  the  simplest  case  to 
deal  with.  However,  any  other  choice  can  be  implemented  while  still  using  the  BW  algorithm.  The  model  is  then  first 
used  to  find  steady-state  solutions  through  a  relaxation  in  time  from  some  arbitrarily  chosen  initial  state  while  holding 
the  boundary  values  of  the  temperature  and  density  constant.  The  steady  state  is  found  with  computational  times 
comparable  to  those  required  to  solve  the  steady-state  equations  alone.  It  is  anticipated  that  any  additions  to  the 
physical  processes  included  in  the  equations  of  motion  would  result  in  the  relaxation  being  the  most  efficient  way, 
computationally,  to  produce  steady-state  solutions. 

The  model  is  next  used  to  give  examples  comparing  the  passage  of  a  transient  on  polytropic  flow  with  the  passage  of 
a  similar  transient  on  thermally  conductive  flow.  This  is  done  for  (i)  initial  states  similar  to  each  other  and  boundary 
conditions  identical  to  those  chosen  by  Whang  and  Chang  (1965)  in  a  very  early  solution  to  the  solar  wind  equations 
with  thermal  conduction,  and  (ii)  similar  initial  states  with  somewhat  lower  base  densities  than  in  the  first  case,  which 
also  results  in  similar  final  states  after  the  base  temperature  has  been  instantaneously  increased  to  introduce  the 
transient.  Large  differences  are  shown  to  exist.  These  differences  are  primarily  due  to  the  ability  of  a  heat  pulse  to 
propagate  at  many  times  the  sound  speed. 

Application  of  the  algorithm  to  different  problems  and  improvements  to  the  model  to  better  describe  the  physics  are 
briefly  outlined  so  as  to  provide  the  necessary  guidance  in  order  to  apply  the  techniques  described  here  to  different 
problems. 
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APPENDIX 

For  the  solar  wind  equations  written  in  matrix  form  (eq.  [15]),  the  vector  of  conserved  variables  V,  the  flux  vectors  F 
and  V,  and  the  vector  of  inhomogeneous  terms  are 


r  i_  2 

~b+rp  . 

^(e  +  rV) 


1=  -2rP  +  b~ ~T 


where  p  is  the  gas  pressure  (eq.  [13]),  T  is  the  temperature  (eq.  [14]),  and  k  is  the  coefficient  of  thermal  conduction  (eqs. 
[6]— [8]).  Using  these  vectors,  the  Jacobian  matrices  A ,  C,  E,  P,  and  P„  defined  in  (19),  are 
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in  which  equations  (A8c,d,e,f)  are  used  in  the  evaluation,  along  with 
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Helium  abundance  (A(Hel)  enhancements  observed  with  Los  Alamos  instruments  on  IMP  6.  7,  and 
8  during  the  1972-1978  interval  have  been  investigated.  Statistical  analysis  of  73  large  events  with 
At  He)  a  10%  (HAEs)  provides  evidence  for  a  close  link  between  helium  enhancements  at  1  A.U.  and 
transient  coronal  mass  ejections.  HAE  events  are  sporadic,  sometimes  clustered  in  time,  and  their 
frequency  of  occurrence  is  approximately  in  phase  with  the  solar  cycle.  Nearly  50%  of  HAEs  are 
associated  with  interplanetary  shocks  and/or  geomagnetic  activity  sudden  commencements,  but  the 
plasma  pattern  associated  with  A(He)  enhancements  is  independent  of  shock  occurrence.  This  pattern 
features  high  magnetic  field  strength,  low  alpha-proton  velocity  difference,  and  low  proton  tempera¬ 
ture.  These  plasma  properties  suggest  that  the  enhancement  is  embedded  in  a  'closed.'  magnetically 
dominated  structure  that  expands  adiabatically.  In  fact,  HAEs  are  likely  to  be  still  evolving 
dynamically  at  I  A.U.  Evidence  of  a  significant  association  between  helium  enhancements  at  I  A.U. 
and  type  II  and  IV  radio  bursts  in  the  corona  is  presented.  We  interpret  these  results  to  mean  that  most 
HAE  events  originate  in  transient  coronal  disturbances  in  which  the  magnetic  field  strength  is 
enhanced.  » 


!.  Introduction 

Helium  abundance  with  respect  to  hydrogen  (hereinafter 
referred  as  A(He))  is  a  parameter  that  varies  greatly  within 
the  solar  wind  at  I  A.U.  Measurements  collected  in  situ  over 
the  last  15  years  [e.g.,  Neugebauer  and  Snyder,  1966; 
Robbins  et  al.,  1970;  Ogilvie  and  Wilkerson,  1969;  Bame, 
1972;  Bollea  et  a!.,  1972;  Feldman  et  al.,  1977;  Neugebauer, 
1981  ]  indicate  that  the  solar  wind  A(He)  varies  from  less  than 
0.001  to  greater  than  0.35.  Such  variability  was  not  expected 
a  priori  and  has  received  considerable  attention  since  its 
discovery,  both  theoretically  [e.g.,  Jokipii,  1966;  Nakada, 
1969;  Joselyn  and  Holzer,  1978;  Geiss  et  al.,  1970;  Holtweg 
etal.,  1978;  Borrini  and  Noci,  1979;  Hollweg,  1981;  Kopp  et 
al..  1981]  and  observationally  [e.g.,  Hirshberg  el  al.,  1972 a, 
b,  1974;  Ogilvie,  1972;  Moreno  and  Palmiotlo,  1973;  Ogilvie 
and  Hirshberg,  1974;  Bame  et  al.,  1977;  Feldman  et  al., 
1977,  1978;  Gosling  et  al.,  1978,  1981;  Borrini  et  al.,  1982], 
Out  of  these  studies  has  come  the  recognition  that  much  of 
solar  wind  A(He)  variability  appears  to  be  related  to  the 
large-scale  structure  of  the  interplanetary  plasma  (i.e.,  high 
speed  streams,  stream  interfaces,  transient  disturbances, 
and  sector  structure). 

In  particular,  flows  with  relative  helium  content  consider¬ 
ably  larger  than  the  average  value  of  5%  have  been  repeated¬ 
ly  associated  with  transient  disturbances.  Some  early  studies 
[e.g..  Gosling  et  al.,  1967;  Bame  et  al.,  1968;  Lazarus  and 
Binsack,  1969;  Ogilvie  et  al.,  1968;  Hirshberg  et  al.,  1971] 
suggested  that  anomalously  high  values  of  A(He)  were  found 
on  occasion  within  a  day  or  so  following  the  passage  of  an 
interplanetary  shock  at  1  A.U.  A  later,  more  thorough, 
statistical  study  of  16  events  where  A(He)  exceeded  0.15 
established  beyond  doubt  that  high  values  of  A(He)  are 
preferentially  (but  not  exclusively)  associated  with  interplan- 
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etary  shocks  and  large  solar  flares  [ Hirshberg  et  al.,  1972a]. 
Further  examples  of  large  A(He)  enhancements  following 
shocks  have  been  presented  in  recent  years  [e.g.,  Bame  et 
al.,  1979,  1981;  Gosling  et  al.,  1980].  In  fact,  about  one  half 
of  all  shocks  observed  at  I  A.U.  are  followed  by  A(He) 
enhancements  >0.8  [Borrini  et  al,  1982].  On  the  other  hand, 
A(He)  enhancements  need  not  be  preceded  by  shocks;  they 
are  often  observed  in  the  low  speed  solar  wind  with  no  shock 
association  [e.g.,  Fenimore,  1980].  Both  with  and  without  a 
shock  association,  A(He)  enhancements  often  have  unusual¬ 
ly  high  ionization  temperatures,  indicative  of  an  origin  in 
active  solar  processes  [ Bame  et  al,  1979;  Fenimore,  1980]. 
On  at  least  one  occasion,  though,  an  A(He>  enhancement 
had  an  unusually  low  ionization  temperature  [ Gosling  et  al, 
1980].  Collectively,  these  observations  suggest  that  A(He) 
enhancements  in  the  solar  wind  signal  the  arrival  of  plasma 
ejected  from  low  in  the  corona  during  a  disturbance  such  as  a 
large  solar  flare  or  an  eruptive  prominence.  It  has  been 
proposed  that  such  plasma  may  be  particularly  rich  in  helium 
because  of  gravitational  settling  in  the  solar  atmosphere 
[Hirshberg  et  al.,  1970]  or  because  of  abnormal  temperature, 
pressure,  or  density  gradients  associated  with  solar  activity 
[Bame  et  al.,  1977]  or  because  of  a  runaway  effect  associated 
with  large  proton  acceleration  [ Borrini  and  Noci.  1979], 
Heretofore,  with  the  exception  of  the  work  by  Fenimore 
[1980],  the  emphasis  of  studies  of  A(He)  enhancements  has 
been  on  shock-related  events.  The  purpose  of  the  present 
work  is  a  thorough  analysis  of  large  helium  enhancements  in 
the  solar  wind,  with  no  restriction  to  shock-related  events. 
Our  study  provides  a  firmer  statistical  basis  for  some  earlier 
results  and  establishes  several  new  associations  for  A(He) 
enhancements  (hereinafter,  following  Fenimore,  referred  to 
as  HAEs).  The  solar  wind  plasma  data  set  used  in  our 
analysis  has  been  collected  by  Los  Alamos  National  Labora¬ 
tory  instruments  about  IMP  6-7-8  from  1971  through  1978. 
One  hour  averages  have  been  employed.  The  instruments 
and  data  reduction  procedures  have  been  discussed  in  previ- 
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ous  publications  [e.g..  As  bridge  et  at.,  1976;  Feldman  et  al., 
1976;  Bame  et  al.,  1977],  and  the  pertinent  details  for  A(He) 
measurements  are  outlined  in  recent  papers  (e.g.,  Borrini  et 
al.,  1981;  Gosling  et  al.,  1981].  Magnetic  field  data  has  been 
obtained  from  the  NSSDC  [King,  1977,  1979]  and  solar  radio 
data  from  the  NOAA  Solar  Geophysical  Data  Book. 

2.  A  Statistical  Analysis  of  Helium  Abundance 
Enhancements 

A(He)  measurements  greater  than  the  long-term  average 
of  5%  occur  frequently,  particularly  in  the  low  speed  solar 
wind  where  A(He)  is  generally  highly  variable  (e.g.,  Bame  et 
al.,  1977;  Feldman  et  al.,  1977].  For  example.  Figure  1  of 
Borrini  et  al.  [1981]  indicates  that  A(He)  is  greater  than  0.07 
approximately  10%  of  the  time.  To  restrict  our  analysis  to  a 
manageable  subset  of  the  data  and  to  eliminate  variations 
that  might  not  be  related  to  the  phenomenon  we  wish  to 
discuss,  we  have  chosen  to  study  HAEs  in  which  A(He) 
>0.10.  Examination  of  the  above-mentioned  figure  reveals 
that  such  events  encompass  —1%  of  all  hourly-averaged 
measurements.  (Note  that  our  choice  is  different  from  that 
used  in  several  previous  studies.  Hirshberg  et  al.  [1972a] 
chose  to  study  events  where  A(He)  >0.15;  Fenimore  [1980] 
studied  events  where  A(He)  >0.09;  and  Borrini  et  al.  (1982] 
discussed  post-shock  events  with  A(He)  20.08.)  To  be 
included  in  the  present  study,  an  event  had  to  have  A(He) 
20.10  for  at  least  2  hours  within  a  24  hour  interval.  The 
starting  time  of  the  HAE  was  chosen  to  be  the  time  of  the 
first  hourly  averaged  value  in  which  A(He)  >0.10.  If  an 
event  persisted  for  many  hours,  but  with  interruptions  in 
which  A(He)  became  less  than  0.10  for  periods  less  than  24 
hours,  the  first  starting  time  was  kept  and  the  event  was 
considered  a  single  HAE.  In  all,  73  HAE  events  were 
identified  in  this  way;  starting  times  for  the  events  were 
tabulated  in  Table  1.  We  are  confident  that  these  73  events 
encompass  the  entire  range  of  large  helium  enhancements  in 
the  solar  wind  and  can  be  used  to  characterize  the  major 
features  of  such  events.  Examination  of  Table  I  reveals  that 
HAEs  are  prevalently  isolated  events,  sometimes  clustered 
in  time,  and  very  rarely  recurrent.  Also,  the  frequency  of 
HAE  observations  is  related  to  the  solar  cycle  (i.e. .  it  is  high 
near  solar  maximum  (1971,  1972,  1977,  1978)  and  low  near 
solar  minimum  (1973-1976)). 

Three  examples  of  our  selected  HAE  events  are  displayed 
in  Figures  1-3.  From  top  to  bottom  these  figures  show  the 
solar  wind  A(He),  proton  density,  bulk  flow  speed,  and 
proton  temperature  for  6  day  intervals  approximately  cen¬ 
tered  on  each  of  three  separate  HAE  events.  One  of  these 
events,  the  HAE  of  May  22-23,  1978  (Figure  3)  was  preced¬ 
ed  by  -I  day  by  an  observed  shock  at  IMP.  while  the  other 
two  events  (December  2,  1974,  and  May  2-4.  1977)  occurred 
independent  of  either  a  shock  or  the  sudden  commencement 
of  a  geomagnetic  storm.  Although  the  HAE  events  of  May  2- 
4,  1977  (Figure  2)  and  May  22-23.  1978  (Figure  3)  have  been 
treated  as  single  events  in  our  study,  they  have  the  appear¬ 
ance  of  multiple  events.  They  may.  in  fact,  have  been 
produced  by  multiple  solar  activity  outbursts.  The  abrupt 
rise  in  A(He)  at  the  onset  of  the  events  in  Figures  1-3  is  a 
characteristic  common  to  many  HAEs.  Finally,  there  is  a 
clear  tendency  for  the  HAE  events  to  occur  at  times  of  low 
proton  temperatures,  a  correlation  consistent  with  previous 
work  [Gosling  et  al.,  1973;  Montgomery  et  al.,  1974;  Bame 
eta!.,  1979], 


Figure  4  presents  the  results  of  a  superposed  epoch 
analysis  of  the  solar  wind  plasma  and  field  data  by  using  as 
key  times  the  onsets  of  our  73  selected  HAE  events.  The 
quantities  plotted  are  A(He),  proton  density,  bulk  speed, 
alpha-proton  speed  difference,  magnetic  field  strength,  and 
proton  temperature.  A  ±6  day  window  about  zero  epoch  has 
been  chosen  so  as  to  allow  an  adequate  assessment  of  the 
statistical  significance  of  the  results.  The  error  bar  plotted  on 
the  middle  of  the  central  line  of  each  plot  is  the  average 
variance  of  the  parameter  for  the  12-day  interval.  The 
average  A(He)  profile  is  asymmetrical,  rising  rapidly  (within 
~4  hours)  to  a  peak  of  0.17  at  zero  epoch,  remaining  above 
0.10  for  several  hours,  and  then  decaying  slowly  back  to 
normal  over  a  period  of  ~  1 .5  days.  This  asymmetry  may  be 
an  artifact  of  the  way  in  which  the  data  were  superimposed; 
however,  the  asymmetry  is  consistent  with  our  examination 
of  individual  events.  In  approximately  half  of  the  HAE 
events  the  initial  rise  in  AfHe)  is  abrupt  and  is  followed  by  a 
gradual  decay  to  average  values  over  1-2  days.  Many  other 
events  consist  of  an  isolated  peak  in  A(He)  lasting  but  a  few 
hours.  In  addition,  some  of  our  HAE  events  are  complex, 
such  as  illustrated  in  Figures  2  and  3,  with  the  secondary 
enhancement  beginning  at  a  variable  lag  relative  to  the  initial 
enhancement. 

Neither  the  proton  density  nor  the  bulk  flow  speed  show 
any  pronounced  average  behavior  in  association  with  helium 
abundance  enhancements.  A  very  small  maximum  in  proton 
density  occurs  —0.5  days  before  the  maximum  in  A(He). 
This  maximum,  which  is  of  marginal  statistical  significance, 
can  be  ascribed  to  interplanetary  compression  preceding  the 
onset  of  helium  enhancements.  For  instance.  44%  of  our 
HAEs  follow  (within  2  days)  the  passage  of  an  interplanetary 
shock  either  observed  directly  or  inferred  from  ground 
magnetograms  (see  Table  I).  On  the  other  hand,  a  density 
spike  of  possible  coronal  origin  often  is  found  at  the  leading 
edge  of  HAE  events  [Bame  et  al.,  1979J. 

A  pronounced  minimum  in  the  average  helium-proton 
speed  difference  is  centered  approximately  on  the  maximum 
in  average  A(He).  Previous  studies  have  indicated  that 
minimums  in  the  speed  difference  are  normally  a  result  of 
enhanced  Coulomb  coupling  between  He  * 1  and  H  ’  in 
regions  of  high  density  and  low  temperature  [ Neugebauer , 
1976],  Although  the  proton  density  is  not  enhanced  at  HAE 
crossing  time,  the  proton  temperature  is  depressed  there, 
suggesting  a  possible  Coulomb  coupling  enhancement.  In 
fact,  using  the  formulas  in  Neugebauer.  one  expects  to 
observe  a  minimum  in  speed  difference  of  the  magnitude 
which  is  observed.  However,  the  minimum  may  be  the  result 
of  other  processes.  For  example,  magnetic  confinement  or 
an  absence  of  wave  deposition,  or  both,  may  account  for  the 
observed  result. 

A  broad,  and  statistically  significant,  maximum  in  magnet¬ 
ic  field  strength  coincides  with  the  maximum  in  A(He). 
However,  the  field  enhancement  begins  —  I  day  prior  to  the 
A(He)  maximum  and  persists  for  —2.5  days  thereafter.  That 
is,  the  HAE  events  on  the  average  occupy  the  core  of 
regions  of  enhanced  magnetic  field  strength.  Although  some 
of  the  field  enhancement  prior  to  the  onset  of  the  HAE  may 
be  the  result  of  compression  in  interplanetary  space,  it  is 
clear  that  the  bulk  of  the  field  enhancement  can  not  be 
explained  in  this  way  because  the  data  show  no  evidence  of 
compressive  heating.  Rather,  a  magnetic  field  enhancement 
is  an  intrinsic  signal  associated  with  HAE  events.  Examples 
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TABLE  I.  Helium  Enhancements  HAEs  1971/1978 


Radio  Bursts 


A(He)  Maxi¬ 

Radio  Bursts 

Type  11  and/or  IV 

mum  Value  Ob¬ 

Interplanetary 

Sudden  Com¬ 

Type  11  and/or  IV 

Observed  Within 

served  Within 

Shock  Associ¬ 

mencement 

Observed  Ahead 

Days  -2  to  -4 

Year 

Month 

Day 

Hour 

the  Event 

ation 

Association 

Within  5  Days 

Prior  to  HAE 

1971 

4 

22 

22 

13.3 

yes 

yes 

yes 

1971 

4 

30 

13 

11.2 

yes 

1971 

5 

17 

18 

24.5 

yes 

yes 

yes 

yes 

1972 

2 

18 

14 

19.8 

yes 

yes 

1972 

2 

26 

17 

12.3 

yes 

yes 

1972 

3 

7 

5* 

13.1 

yes 

yes 

yes 

yes 

1972 

3 

9 

16* 

23.4 

yes 

yes 

1972 

3 

II 

8* 

11.5 

yes 

yes 

1972 

3 

27 

21 

13.2 

1972 

4 

13 

22* 

10.4 

yes 

yes 

1972 

4 

16 

4* 

10.7 

1972 

4 

17 

20* 

11.5 

yes 

1972 

4 

20 

7* 

12.5 

1972 

4 

21 

13* 

111 

yes 

yes 

1972 

4 

26 

12* 

13.1 

1972 

4 

29 

5* 

19.9 

1972 

5 

16 

2 

11.2 

yes 

yes 

yes 

yes 

1972 

6 

21 

20 

14.5 

yes 

1972 

II 

29 

8 

33.9 

yes 

yes 

1973 

1 

19 

1 

11.0 

1973 

1 

20 

18 

14.8 

1973 

2 

18 

18 

11.8 

yes 

yes 

1973 

4 

13 

17 

16.1 

yes 

yes 

yes 

yes 

1973 

5 

3 

8* 

15.1 

yes 

yes 

yes 

1973 

5 

5 

1* 

17.8 

yes 

yes 

1973 

5 

21 

13 

12.6 

yes 

yes 

yes 

yes 

1973 

5 

25 

6 

14.2 

yes 

yes 

1973 

9 

13 

7 

11.6 

yes 

yes 

1973 

9 

26 

16 

12.7 

yes 

yes 

1974 

1 

25 

0 

10.8 

1974 

2 

22 

12 

10.8 

yes 

yes 

1974 

5 

8 

16 

20.5 

yes 

1974 

6 

10 

22 

11.5 

yes 

yes 

yes 

yes 

1974 

7 

6 

14 

20.5 

yes 

yes 

yes 

yes 

1974 

7 

31 

16 

11.7 

yes 

1974 

9 

16 

14 

14.5 

yes 

yes 

yes 

yes 

1974 

10 

13 

15* 

14.5 

yes 

yes 

yes 

yes 

1974 

10 

15 

23* 

12.0 

yes 

yes 

yes 

1974 

12 

2 

12 

14.3 

yes 

yes 

1975 

1 

6 

II 

12.0 

1975 

3 

7 

16 

22.4 

yes 

1976 

5 

3 

8 

13.3 

yes 

yes 

yes 

1976 

12 

18 

21 

21.25 

1977 

2 

14 

16 

16.0 

yes 

yes 

1977 

5 

2 

21 

19.5 

1977 

5 

29 

16 

1 1.0 

yes 

1977 

7 

5 

7 

39.4 

1977 

9 

20 

10* 

11.3 

yes 

yes 

1977 

9 

22 

19* 

12.5 

yes 

yes 

yes 

yes 

1977 

9 

27 

7 

10.7 

yes 

yes 

1977 

10 

4 

18 

10.3 

yes 

1977 

10 

21 

2 

40.2 

1977 

10 

27 

23 

21.6 

yes 

yes 

1977 

II 

5 

6* 

15.7 

1977 

II 

6 

8* 

III 

1977 

12 

15 

6 

12.4 

yes 

1978 

1 

5 

17 

19.9 

yes 

yes 

yes 

yes 

1978 

1 

30 

13 

36.1 

yes 

1978 

3 

5 

6 

13.9 

1978 

4 

10 

23 

29.2 

yes 

yes 

yes 

yes 

1978 

5 

4 

12 

22.9 

yes 

yes 

1978 

5 

8 

II 

16.3 

yes 

yes 

1978 

5 

22 

18 

21.5 

yes 

yes 

yes 

yes 

1978 

5 

31 

15 

18.2 

yes 

yes 

1978 

6 

4 

9 

1 1.0 

yes 

yes 

yes 

yes 

1978 

7 

5 

18 

16.3 

yes 

yes 

yes 

yes 

1978 

7 

16 

0* 

21.2 

yes 

yes 

yes 

1978 

7 

18 

II* 

14.1 

yes 

yes 

1978 

9 

8 

15* 

10.0 

yes 

yes 

yes 
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TABLE  I  ( continued l 


Year 

Month 

Day 

Hour 

A(He)  Maxi¬ 
mum  Valife  Ob¬ 
served  Within 
the  Event 

Interplanetary 
Shock  Associ¬ 
ation 

Sudden  Com¬ 
mencement 
Association 

Radio  Bursts 
Type  II  and/or  IV 
Observed  Ahead 
Within  5  Days 

Radio  Bursts 
Type  II  and/or  IV 
Observed  Within 
Days  -2  to  -4 
Prior  to  HAE 

1978 

9 

10 

9* 

12.7 

yes 

yes 

yes 

1978 

9 

30 

5 

37.4 

yes 

yes 

yes 

yes 

1979 

1 

1 

16* 

13.5 

yes 

yes 

1979 

1 

3 

13* 

11.2 

yes 

yes 

‘Possible  multiple  event. 

Note:  The  association  with  interplanetary  shocks  or  sudden  commencements  refer  to  4K  hours  prior  to  the  event. 


of  noncompressive  magnetic  field  enhancements  have  also 
been  reported  by  Burlaga  and  Kina  [  1979|  and  Klein  and 
Burlaga  [1982). 

Finally,  a  deep  minimum  in  average  proton  temperature 
coincides  with  the  HAE.  This  minimum  is  expected  from  our 
discussion  of  Figures  1-3  and  from  previous  work.  Although 
the  temperature  depression  and  the  HAE  overlap  in  Figure 
4,  an  examination  of  individual  events  reveals  occasions 
when  they  are  spatially  distinct.  Because  HAE  events  often 
exhibit  anomalously  hot  ionization  temperatures  [Bame  et 
ai.  1979;  Fenimore,  1980),  it  is  unlikely  that  the  low  proton 
temperature  is  of  coronal  origin.  A  more  likely  explanation  is 
that  these  low  temperatures  result  from  enhanced  cooling  in 
interplanetary  space  [e.g..  Gosling  et  ai,  1973;  Montgomery 


TIME  (days)  DEC  1974 

Fig.  I.  One  hour  averaged  solar  wind  data.  November  30  to 
December  4.  1974.  Notice  the  noncompressive  density  enhancement 
preceding  the  HAE. 


et  al..  1974)  produced  by  either  an  adiabatic  expansion 
within  a  closed  magnetic  field  topology  or  a  faster  than 
normal  expansion  of  the  HAE  plasma,  or  both. 

Figure  5  shows  the  result  of  a  superimposed  epoch  analy¬ 
sis  of  total  pressure  IP/  =  i.n.AT,  +  B2/Kn).  gas  pressure 
(±,n,k  T, )  and  magnetic  pressure  (B:iXrr)  where  /'  =  1,3  for 
electrons,  protons,  and  alpha  particles,  respectively;  k  is 
Boltzman's  constant;  and  B  is  the  magnetic  field  strength) 
keying  as  in  Figure  4  on  onsets  of  our  HAE  events.  A  broad 
maximum  in  total  pressure  surrounds  the  onset  time,  al¬ 
though  the  maximum  in  A(He)  coincides  with  what  appears 
to  be  a  local  minimum  in  total  pressure.  The  double  peaked 
nature  of  the  pressure  profile  results  from  having  a  high 
proton  density  and  high  field  strength  prior  to  the  HAE.  a 


Fig.  2.  One  hour  averaged  solar  wind  data.  April  30  to  May  5. 
1977.  Notice  the  multi-peak  structure  of  the  helium  enhancement 
and  the  contemporary  proton  temperature  depression. 
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Fig.  3.  One  hour  averaged  solar  wind  data.  May  20-25,  1978. 
An  interplanetary  shock  was  observed  on  May  21 .  Notice  the  multi¬ 
peak  structure  of  the  helium  enhancement. 


minimum  in  temperature  coincident  with  the  HAE.  and  an 
enhanced  field  strength  in  the  decaying  phase  of  the  HAE. 
The  overall  pressure  maximum  associated  with  HAE  events 
implies  that,  on  the  average.  HAEs  must  be  part  of  strut 
tures  still  dynamically  evolving  in  interplanetary  space. 

Examination  of  Table  I  reveals  that  33*/?  of  our  HAE 
events  are  preceded  or  followed  within  ±3  days  by  a 
companion  HAE.  In  view  of  the  fact  that  only  73  events  have 
been  identified  in  nearly  8  years  of  relatively  continuous  data 
coverage,  this  indicates  that  HAE  events  tend  to  occur  in 
clusters.  Furthermore,  because  of  our  definition  of  HAE 
events,  the  above  percentage  does  not  fully  represent  the 
phenomenon.  Many  examples  of  complex  events  such  as 
displayed  in  Figures  2  and  3  are  present  in  the  data.  Although 
such  events  are  included  in  our  study  as  single,  they  may  in 
fact  be  multiple  events.  (On  the  other  hand,  this  structuring 
may  also  be  intrinsic  to  single  HAE  events  [Bame  el  al., 
1979)).  A  similar  clustering  tendency  for  interplanetary 
shocks  has  been  reported  previously  [ Borrini  et  al.,  1982). 
We  believe  that  both  multiple  shocks  and  multiple  HAE 
events  are  caused  by  multiple  outbursts  from  the  same  solar 
activity  center  on  a  time  scale  of  several  days  as  reported, 
for  example,  by  Hildner  et  al.  1 1976). 

It  is  not  obvious  from  an  examination  of  individual  HAEs 
that  there  is  anything  fundamentally  different  in  the  charac¬ 
ter  of  events  associated  with  shock  wave  disturbances  as 
opposed  to  those  not  so  associated.  To  test  this  hypothesis, 
we  repeated  our  superposed  epoch  analysis  separating  the 
events  according  to  their  association  or  lack  of  association 
with  shock  events  (and  sudden  commencements).  The  only 
appreciable  difference  in  the  resulting  profiles  is  that  the  flow 
speed  is  somewhat  higher,  on  the  average,  at  the  onset  of  the 
HAE  for  the  shock-associated  events. 


KEYS  =  73  HELIUM  ENHANCEMENT  EVENTS 

A(Hv)  >10%  FOR  MORE  THAN  2  HOURS  WITHIN  A  DAY 


Fig.  4.  Superposed  epoch  plots  of  solar  wind  alpha-proton  flux  ratio  (helium  abundance),  proton  density,  bulk 
speed,  alpha-proton  velocity  difference,  magnetic  field  intensity,  and  proton  temperature  for  73  helium  abundance 
enhancements  at  I  A.U.  (For  lack  of  data  only  42  keytimes  are  used  in  the  magnetic  field  strength  plot.)  One  hour 
averages  of  data  are  employed  and  the  first  hour  with  A(He)  2|09f  is  defined  as  the  starting  time  of  the  enhancement. 
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Fig.  5.  Superposed  epoch  plots  of  solar  wind  total  pressure,  gas 
pressure,  and  magnetic  pressure  for  the  same  key  time  set  of  7.3 
helium  abundance  enhancements  of  Figure  4.  Again  because  of 
missing  data,  only  42  key  times  are  used  in  the  magnetic  pressure 
plot,  69  in  the  gas  pressure  plot,  and  40  in  the  total  pressure  plot. 


3.  The  Association  of  Helium  Abundance 
Enhancements  With  Type  11 
and  IV  Solar  Radio  Bursts 

As  was  already  mentioned,  approximately  half  of  our 
HAE  events  follow  within  2  days  of  a  transient  interplane¬ 
tary  shock  wave  disturbance.  Moreover,  because  there 
appears  to  be  nothing  fundamentally  different  between  shock 
and  non-shock  associated  HAE  events,  it  is  the  authors' 
opinion  that  most,  if  not  all.  HAE  events  detected  at  I  A.U. 
can  be  associated  with  the  kind  of  transient  mass  ejection 
events  so  spectacularly  photographed  at  the  solar  limbs  by 
orbiting  coronagraphs  [e.g..  Stewart  et  al..  1974:  MacQueen 
et  al..  1974],  However,  few  HAE  events  observed  by  earth- 
orbiting  satellites  would  be  expected  to  show  a  direct 
association  with  the  limb  events  photographed  by  corona- 
graphs.  To  test  for  a  definitive  association  of  HAEs  with 
coronal  mass  ejection  events,  we  have  studied  the  connec¬ 
tion  between  HAE  events  and  decametric  type  II  and  type 
IV  solar  radio  bursts.  Such  bursts  have  been  extensively 
studied  by  solar  radio  astronomers  [see.  for  example.  Smerd 
and  Dulk,  1971;  Wild  and  Sttierd.  1972;  Stone,  1974).  are 
strongly  correlated  with  large  solar  flares  and  eruptive 
prominences,  and  are  believed  to  be  manifestations  of  shock 
disturbances  in  the  corona.  Type  II  and  IV  bursts  also  have  a 
very  strong  correlation  with  coronal  mass  ejection  events: 
virtually  all  such  bursts  originating  near  the  solar  limb  during 
the  Skylab  mission  were  accompanied  by  coronal  mass 
ejections  observed  by  the  coronagraph  [Goslinn  et  al..  1974: 
Munro  et  al..  I979|.  However,  only  30 c/r  of  all  mass  ejection 
events  observed  during  this  same  period  produced  type  II 
and/or  IV  radio  bursts.  In  particular,  it  was  the  fastest 
ejections  (observed  speed  -400  km)  that  produced  such 
bursts  \Go\lmk.  et  al..  I976|.  In  summary,  type  II  and/or  IV 


radio  bursts  are  excellent  indirect  evidence  of  fast  mass 
ejection  events  in  the  solar  corona. 

Using  the  reports  of  decametric  type  II  and  IV  bursts  of 
the  Solar  Geophysical  Data  Book  distributed  by  NOAA.  one 
finds  that  68%  of  our  HAE  events  are  preceded  wit'-m  >  days 
by  a  reported  type  II  and/or  type  IV  radio  burst;  78%  of  the 
shock-associated  HAE  events  fit  in  this  category  as  well  as 
61%  of  the  nonshock-associated  events.  If  we  restrict  the 
interval  of  correlation  to  2-4  days  prior  to  the  HAE  event 
onset,  59%  of  our  events  are  preceded  by  type  II  and/or  IV 
radio  bursts  (72%  of  the  shock-associated  and  40%  of  the 
nonshock-associated  events).  These  percentages,  however, 
do  not  adequately  reflect  the  degree  of  correlation  between 
HAE  events  and  type  II  and  IV  radio  bursts  because  they 
say  nothing  about  the  average  occurrence  frequency  of  the 
radio  bursts. 

Figure  6  shows  a  plot  of  the  total  number  of  type  II  and  IV 
radio  bursts  reported  as  a  function  of  time  relative  to  the 
observation  of  our  73  HAE  events  at  1  A.U.  (day  zero).  The 
plot  includes  at  least  4  lag  days  (days  -12,  -7.  +1,  tt| 
sufficiently  removed  from  the  time  of  observation  of  the 
HAE  events  (day  zero)  as  to  allow  an  estimate  of  the  random 
occurrence  of  type  II  and  IV  bursts  in  the  corona.  The 
number  of  random  occurrences  averaged  over  the  above  4 
lag  days  is  52,  corresponding  to  a  random  occurrence  rale  of 
0.71  per  day.  On  the  other  hand,  the  average  number  of 
bursts  reported  2-4  days  prior  to  onset  of  the  HAE  events 
(where  the  distribution  of  bursts  peaks)  is  128.  correspond¬ 
ing  to  an  occurrence  rate  2.5  times  greater  than  random.  A 
lag  of  2-4  days  is  consistent  with  a  reasonable  range  of 
average  sun-earth  transit  speeds  (435-870  km  s  Howev¬ 
er,  not  only  is  the  distribution  of  bursts  peaked  at  lags  of  2-4 
days,  it  is  also  greater  than  random  at  all  lags  Irom  0  to  6 
days.  This  effect  can  almost  certainly  be  attributed  to  the 
fact  that  radio  bursts,  like  coronal  mass  ejection  events, 
interplanetary  shock  disturbances,  and  HAE  events,  tend  to 
cluster  within  intervals  of  several  days. 

Figures  7  and  8  show  identical  analyses  for  type  IV  and  II 


Fig.  6.  A  histogram  of  the  total  number  of  type  II  and  IV 
decametric  radio  bursts  observed  on  the  sun  as  a  function  of  time 
relative  to  the  onset  of  a  HAE  event  at  I  A.U.  The  peak  of  the 
distribution  coincides  with  the  time  lag  of  2-4  days  (corresponding 
to  an  average  sun-earth  transit  speed  between  435  and  870  km/s). 
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Fig.  7.  Same  as  Figure  ft  for  the  total  number  of  type  IV  radio 
bursts  only. 


bursts  considered  separately.  The  same  correlation  is  pres¬ 
ent  for  each  individual  set  as  for  the  combined  set  (Figure  6), 
but  is  considerably  better  for  type  IV  bursts  than  for  type  II 
bursts,  the  peak  signal  to  noise  ratios  being  about  4.0  and 
1.4,  respectively.  It  is  not  presently  clear  why  the  correlation 
is  better  for  type  IV  bursts  since  presumably  both  types  of 
bursts  are  manifestations  of  fast  outward-moving  coronal 
mass  ejection  events. 

The  foregoing  documents  that  a  correlation  exists  between 
some  HAE  events  and  type  II  and  IV  solar  radio  bursts,  and 
by  inference  fast-moving  coronal  mass  ejection  events. 
However,  before  leaving  this  topic,  a  few  qualifications  are 
in  order.  First,  a  fair  fraction  of  the  radio  burst  events 
included  in  our  statistics  and  in  Figures  6-8  are  almost 
certainly  limb  events,  which  probably  have  little,  if  any, 
direct  effect  upon  solar  wind  conditions  near  the  earth.  Thus, 
for  example,  our  statement  that  59%  of  all  HAE  events  are 
preceded  2-4  days  earlier  by  one  or  more  type  II  and/or  type 
IV  bursts  almost  certainly  overestimates  the  percentage  of 
HAE  events  which  actually  are  directly  associated  with 
coronal  transients  that  produced  such  radio  bursts.  Second, 
we  have  not  established  a  one-to-one  correlation  between 
HAE  events  and  either  type  II  and/or  IV  radio  bursts  or  fast 
coronal  mass  ejection  events.  Nor  is  there  any  good  reason 
to  expect  that  a  one-to-one  correlation  should  in  fact  exist. 
Rather,  we  expect  that  HAE  events  in  the  solar  wind  at  I 
A.U.  are  a  signature  (and  not  necessarily  the  only  one  or  an 
invariant  one)  of  coronal  mass  ejection  events  pretty  much 
independent  of  the  speed  of  the  ejection  event  (see  also 
Ogilvie  and  Hirshberg.  1974).  The  very  fast  ejection  events 
that  produce  type  II  and/or  IV  radio  bursts  are  a  special 
subset  of  such  events.  It  is  important  to  document,  as  has 
been  done  here,  that  this  special  subset  of  mass  ejection 
events  has  a  strong,  positive  correlation  with  HAE  events  at 
I  A.U.  because  it  makes  plausible  the  inference  that  all  HAE 
events  are  caused  by  transient  processes  in  the  solar  atmo¬ 
sphere. 

4.  Observational  Summary 

In  the  previous  sections  we  presented  several  observa¬ 
tions  supporting  the  association  between  HAEs  at  I  A.U. 
and  transient  mass  ejections  in  the  corona.  These  observa¬ 
tions  are  summarized  briefly  below: 


1.  Helium  enhancements  are  sporadic  phenomena, 
sometimes  clustered  in  time,  very  rarely  recurrent. 

2.  The  occurrence  frequency  of  HAEs  is  related  to  the 
solar  cycle  (frequent  at  maximum,  scarce  at  minimum). 

3.  The  HAE  events  defined  in  section  2  are  associated 
with  interplanetary  shocks  in  44%  of  the  cases. 

4.  The  plasma  pattern  associated  with  HAEs  is  indepen¬ 
dent  of  shock  occurrence. 

5.  This  plasma  pattern  features  high  magnetic  field 
strength,  low  proton  temperature,  low  alpha-proton  velocity 
difference,  and  high  total  pressure. 

6.  Helium  enhancements  at  I  A.U.  are  often  associated 
with  type  II  and  IV  radio  bursts  in  the  corona. 

5.  Discussion 

Large  enhancements  in  the  relative  abundance  of  helium 
in  the  solar  wind  occur  relatively  infrequently.  This  paper 
has  concentrated  on  enhancements  wherein  A(He)  achieved 
values  greater  than  0.10,  accounting  for  ~1%  of  all  hourly 
averaged  measurements  in  the  1971-1978  interval.  Increases 
in  A(He)  of  this  magnitude  were  identified  in  73  separate 
events,  only  about  half  of  which  were  preceded  within  2  days 
by  an  observed  shock  or  sudden  commencement  of  a  geo¬ 
magnetic  storm.  Said  another  way,  a  substantial  fraction  of 
solar  wind  helium  abundance  enhancements  occur  in  the 
absence  of  any  shock  wave  disturbance.  These  enhance¬ 
ments,  though,  do  not  seem  to  differ  in  any  substantial  way 
from  the  shock  related  ones;  in  particular  both  types  of 
enhancements  are  associated  with  very  similar  plasma  char¬ 
acteristics. 

A  superposed  epoch  analysis  of  plasma  and  field  data  has 
been  performed  to  determine  the  average  character  of  A(He) 
enhanced  plasma.  Typically,  these  enhancements  occur 
within  plasma  of  average  particle  density  and  flow  speed,  but 
unusually  low  proton  temperature  and  high  magnetic  field 
strength.  Whereas  on  the  average  the  abundance  enhance¬ 
ment  coincides  spatially  with  the  region  of  anomalously  low 
temperature,  the  region  of  high  field  strength  is  considerably 
broader  than  the  abundance  enhancement.  There  is  little 
evidence  of  a  density  compression  or  heating,  which  might 
suggest  that  the  field  enhancement  is  produced  by  dynamical 
processes  in  interplanetary  space.  We  thus  interpret  the  high 
field  measurement  to  mean  that  helium  abundance  enhance¬ 
ments  originate  in  coronal  regions  of  unusually  high  field 
strength.  On  the  other  hand,  we  interpret  the  low  tempera¬ 
tures  to  be  evidence  that  these  magnetic  structures  are  either 
topologically  closed  (hence  adiabatically  cooling)  or  have 
expanded  considerably  faster  than  normal  in  transit  from  the 
sun.  Further  possible  evidence  for  a  closed  structure  sur- 


Fig.  8,  Same  as  Figure  ft  for  the  total  number  of  type  If  radio 
bursts  only. 
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rounding  the  helium  rich  plasma  is  the  observation  of  low 
average  alpha-proton  velocity  difference  within  HAEs.  A 
magnetic  confinement  could  help  to  create  a  region  in  which 
the  process  responsible  for  the  preferential  acceleration  of 
alpha  particles  is  not  effective,  although  this  explanation  is 
not  necessary  [see,  e.g.,  Neugebauer,  1976). 

We  have  presented  evidence  that  suggests  (but  does  not 
prove)  that  HAEs  are  a  1  A.U.  signal  of  coronal  transients. 
Because  the  typical  coronal  mass  ejection  event  has  a  speed 
at  four  solar  radii  (—470  kms'1)  [Gosling  et  al.,  1976],  which 
is  not  far  different  from  the  average  solar  wind  speed  at  1 
A.U.  (—450  km  s_l),  we  should  not  be  surprised  that  helium 
abundance  enhancements  typically  occur  within  plasma  of 
only  average  speed.  By  way  of  contrast,  one  of  the  charac¬ 
teristics  that  distinguishes  coronal  mass  ejection  events  in 
the  corona  is  their  higher  than  average  density.  The  fact  that 
these  events  seem  to  have  only  average  density  at  I  A.U. 
suggests  that  they  expand  in  transit  from  the  sun  to  a  greater 
degree  than  does  the  ordinary  solar  wind  stream  tube.  This 
suggestion  is  consistent  with  the  observation  that  abundance 
enhancements  typically  are  also  regions  of  high  pressure  (but 
low  proton  temperature)  still  evolving  dynamically  even  at  I 
A.U. 

It  is  unlikely  that  helium  abundance  enhancements  where¬ 
in  A(He)  sO.IO  account  for  all  coronal  mass  ejection  events 
in  the  solar  wind  at  I  A.U.,  since  these  enhancements  occur 
much  less  frequently  than  do  coronal  mass-ejection  events 
that  occur  somewhere  on  the  sun  at  a  rate  of  about  I  per  day. 
depending  on  the  phase  of  the  solar  cycle  [Hildner  et  al., 

1976] .  What,  then,  is  the  I  A.U.  signature  of  these  other 
mass  ejection  events,  and  why  should  the  signature  vary? 
One  possibility  is  that  a  subset  of  abundance  enhancements 
with  0.05  s  A(He)  s  0. 10  also  signal  mass  ejection  events  at 
1  A.U.  Another  possibility  is  that  a  subset  of  so-called 
'noncompressive  density  enhancements’  [Gosling  et  id 

1977]  are  such  a  signal  [e.g.,  Fenimore,  1980|.  (On  the  other 
hand,  a  large  fraction  of  noncompressive  density  enhance¬ 
ments  have  anomalously  low  A(He)  and  almost  certainly 
signal  spacecraft  intersections  with  coronal  streamers  [Bor¬ 
rini  et  al..  1981:  Gosling  et  al.,  I98l|.  And,  finally,  it  is 
possible  that  magnetic  ‘clouds’  [e.g.,  Burlaga  and  Behan- 
non,  I982|  without  A(He)  enhancements  form  a  subset  of 
coronal  mass  ejection  events  at  I  A.U.  The  great  variety  of 
physical  appearances  of  mass  ejection  events  in  the  corona 
[see.  for  example,  the  discussion  by  Munro  et  al..  1979] 
suggests  that  the  physical  structure  of  these  events  may  vary 
considerably  from  one  event  to  the  next.  Altogether,  then, 
there  is  no  reason  to  expect  that  a  single  anomaly  such  as  an 
A(He)  enhancement  would  be  sufficient  to  identify  all  cor¬ 
onal  mass  ejections  in  the  solar  wind.  On  the  other  hand,  on 
the  basis  of  previous  and  our  present  findings,  we  feel 
confident  that  major  helium  enhancements  are  one  of  the 
clearest  indications  of  transient  mass  ejections  at  I  A.U. 
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